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NOTATION 


ROMAN 


1. Italicised Lower Case 


a 

b 

d 

f 

g 

k 

1 

m 

<7 

r 

t 

v 

x 

y 


constant 

constant 

diameter 

frequency 

scalar mass flux 

turbulence kinetic energy 

length 

index limit 

contact heat flux scalar 

radius 

time 

velocity scalar 

displacement scalar or variable 
variable 


2. Italicised Upper Case 


A 

C 

E 

G 

H 

K 

M 

N 

P 

R 

T 

U 

V 


area 

heat capacity 

external and mutual energy 

total flux 

enthalpy 

cons tant 

mass 

non-dimensional parameter 

pressure 

gas constant 

temperature 

internal energy 

volume 


ill 



3. Bold Lower Case 


f external and mutual force 

g mass flux density 

n unit outward normal 

q contact heat flux 

s contact force 

v velocity 

x displacement 


4, Bold Upper Case 

D deformation tensor 

I identity matrix 

T extra stress tensor 


5. Computer Programme Variables 


NHL 

NKF 

NRF 

NRL 


index of 
index of 
index of 
index of 


last heater mass/energy discrete volume 
first cooler mass/energy discrete volume 
first regenerator mass/energy discrete volume 
last regenerator mass/energy discrete volume 


GREEK 

1 , Upper Case 

A bulk compressibility 

$ dissipation 


2 . Lower Case 


a constant 

P constant 

S unit vector component 

€ turbulence kinetic energy dissipation rate 


1 V 


p 

V 

p 

T 

(j) 


dynamic viscosity 
kinematic viscosity 
density 

extra shear stress tensor component 
generalised scalar, vector, or tensor quantity 
specific dissipation rate 


HEBREW 

ft transformation tensor 


OPERATORS 


d 

ddt 

D 

f() 

30 

a 

A 

v 

/ 

Z 


total derivative 

total derivative with respect to time 

substantive derivative 

function of 

quantisation function 

partial derivative 

incremental change 

divergence 

integral 

summation 


time average of 

[V j \j> volume average of \j> 


time average of volume average of ^ 
i/> — [tt]^ time average of time average of }f> 


1*1 


n 

u 

c 


absolute value or magnitude of 

scalar product of vectors, vector product of vector and 
tensor 

scalar product of tensors 

intersection 

union 

proper subset 


V 



SUBSCRIPTS 


a 

ch 

i 

I 

j 

k 

M 

(m) 

n 

N 

nL 

nr 

nR 

nx 

P 

Pe 

Pr 

r 

Re 

(s) 

T 

Va 

x 

At 


acoustic 

characteristic 

index 

index 

index 

index 

index limit in two-dimensional space 

material body 

momentum discrete volume 

index limit in two-dimensional space 

momentum discrete volume, left hand 

radial momentum discrete volume 

momentum discrete volume, right hand 

axial momentum discrete volume 

at constant pressure 

Peclet 

Prandtl 

regenerator 

Reynolds 

system of particles 
at constant temperature 
Valensi 
axial 

time increment 


SUPERSCRIPTS 

(t) turbulent 

T transpose 

s previous time step 

* distinguishing indicator 

' fluctuating component 

per unit mass 

• time rate of change 


Vi 



CHAPTER 1 


INTRODUCTION 


The activities described in this report do not constitute a continuum 
but rather a series of linked smaller investigations in the general area of 
one- and two-dimensional Stirling machine simulation. The initial impetus for 
these investigations was the development and construction of the Mechanical 
Engineering Test Rig (METR) under a grant awarded by NASA to Dr Terry Simon at 
the Department of Mechanical Engineering, University of Minnesota. The 
purpose of the METR is to provide experimental data on oscillating turbulent 
flows in Stirling machine working fluid flow path components (heater, cooler, 
regenerator, etc.) with particular emphasis on laminar/turbulent flow 
transitions . 

Hence, the initial goals for the grant awarded by NASA were, broadly, 
to provide computer simulation backup for the design of the METR and to 
analyze the results produced. This was envisaged in two phases: first, to 
apply an existing one -dimensional Stirling machine simulation code to the METR 
and second, to adapt a two-dimensional fluid mechanics code which had been 
developed for simulating high Rayleigh number buoyant cavity flows to the 
METR. The key aspect of this latter component was the development of an 
appropriate turbulence model suitable for generalised application to Stirling 
simulation. A final step was then to apply the two-dimensional code to an 
existing Stirling machine for which adequate experimental data exist. 

The work described herein was carried out over a period of three years 
on a part-time basis. Forty percent of the first year's funding was provided 
as a match to the NASA funds by the Underground Space Center, University of 
Minnesota, which also made its computing facilities available to the project 
at no charge. 


1.1 OBJECTIVES 


With the advantage of a posteriori clarity, the following overall 
objectives guided the course of the work: 

1. Apply an existing one -dimensional simulation code to the METR. 

2. Adapt and apply an existing two-dimensional fluid mechanics code 
to the METR. 

3. Use the METR experimental results to guide the development of a 
turbulence model appropriate for generalised application to 
Stirling machine simulation. 
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4 . 


Validate the two-dimensional simulation including its turbulence 
model against experimental data for an existing Stirling engine. 


1.2 PROJECT EVOLUTION 


The basic simulation codes applied to the METR are the fully implicit, 
discrete volume simulations developed during the course of the author's PhD 
research program (Go87). The application of the one -dimensional version of 
the code to the METR was elementary and provided some design guidance to Simon 
and Seume (the graduate student conducting the METR research) in developing 
the final design of the rig. 

Because a significant delay prior to the commissioning of the rig was 
expected, it was decided to proceed with applying the one- and two-dimensional 
simulations to a Stirling engine, the latter simulation initially without the 
inclusion of a turbulence model. NASA chose the Space Power Demonstrator 
Engine (SPDE) as the target engine. This back-to-back, free-piston design is 
characterised by an operating frequency of 100 Hz and a mean pressurisation of 
150 bars which, combined with a relatively short working fluid flow path, 
yield an engine characteristic number (N c ^) (see section 2.6) of about 25. 

This means that there are roughly 25 complete information propagation 
traverses between the expansion and compression space pistons during each 
cycle. This may be compared with a typical characteristic number of 96 for 
the GM-GPU3 kinematic engine. In the case of the GM-GPU3 engine, the 
characteristic number proved to be large enough so that modelling of 
information propagation effects did not prove necessary in order to match the 
experimental data (Go87) . However, In the case of the SPDE, such modelling of 
information propagation did enable agreement between the measured and 
simulated piston indicated works to be obtained. 

This elicited some controversy not only within NASA but also among 
other Stirling engine analysts who expressed doubt whether information 
propagation effects are physically relevant at low Mach numbers. Hence, a 
significant deviation of the grant was initiated into an investigation of the 
information propagation phenomenon in the context of a physical application 
remote from Stirling engines. While not settling the controversy, this 
investigation did reveal the limitations of the simulation code by 
establishing a lower characteristic number limit below which the code was 
judged inapplicable. 

Within these limits, the two-dimensional fluid dynamics code was 
modified and successfully applied to the SPDE. The methodology adopted was to 
treat the heater as a two-dimensional entity represented by a single "typical" 
tube in an otherwise one -dimensional system. This was intended to facilitate 
the application of the METR turbulence data to an actual engine (via the 
simulation code) since the rig would be configured to represent just such an 
SPDE heater tube. A significant aspect of this application was the 
development of a mesh generation scheme enabling a "seamless" junction between 
the one -dimensional rectilinear and two-dimensional cylindrical spatial 
discretisations . 
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At this stage, the METR still did not yield quantitative data suitable 
for turbulence model development. Nevertheless, the two-dimensional SPDE code 
was ported to the METR geometry, albeit without the inclusion of a turbulence 
model. Some mean velocity qualification comparisons were made between the 
simulated and preliminary experimental data. These comparisons revealed 
systematic errors in the experimental data. The errors included flow 
anomalies (apparently caused by piston/cylinder "sticking"), the absence of 
ambient boundary conditions and hot-wire anemometer calibration uncertainties. 
Qualitatively useful data did however become available when Joerge Seume 
published his PhD thesis (Se88) . Hence, since the grant period was drawing to 
a close, development of the turbulence model was initiated without the benefit 
of an experimental benchmark. Finally, a single half -cycle of quantitatively 
useful turbulence data (despite the continued existence of systematic errors) 
for the SPDE heater tube configuration was eventually delivered about three 
weeks prior to the termination of the grant. This allowed barely enough time 
to make some preliminary turbulence model evaluations and to define the 
critical issues in oscillating flow turbulence modelling. The larger 
objective of applying the turbulence model to the SPDE could not be fulfilled. 

The structure of this report thus reflects the modus operandi of the 
grant itself, comprising a compendium of sectional reports delivered to NASA 
at the termination of each phase. Chapter 2 summarizes the simulation model 
and its theoretical foundations; chapter 3 discusses the simulation of the 
SPDE; chapter 4 describes an investigation of the information propagation 
issue based upon an analytic description of a transmission line; and chapter 5 
is devoted to the simulation of the METR and the development of a turbulence 
model. Salient conclusions and some directions for future research arising 
therefrom are presented chapter 6. 
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CHAPTER 2 


THE SIMULATION MODEL 


2.1 INTRODUCTION 


This chapter presents a narrative overview of the physics and numerics 
of the simulation model. Included is the postulational basis from which the 
conservation equations are derived as well as the definitions of the 
discretisation, turbulence and information propagation models. Since the 
symbolic development of the simulation model is voluminous, the reader is 
directed toward reference Go87 for a complete and rigorous derivation of all 
the equations presented. 

The philosophical basis upon which the simulation model rests is 
described by Tisza (Ti66) as the 'postulational' approach in his discussion of 
the evolution of the concepts of thermodynamics. In summarizing the efficacy 
of the postulational approach, Tisza makes the following critical observation: 

'First, and most important, we claim no absolute validity for our 
postulational basis. The validity of the postulates and the usefulness 
of the primitive concepts are only tentative and have to be justified 
by the experimental verification of the implications of the theory. ' 

Thus the postulational approach used to develop a symbolic description 
of the fluid dynamics of Stirling cycle machines ultimately can be justified 
only by the extent to which the results produced can be given validity by 
experimental observation. 


2.2 THE INTEGRAL CONSERVATION BALANCES 


The integral conservation balances forming the backbone of the 
simulation codes are derived from four postulates. The first postulate is 
based on the classical concept that matter is uniformly distributed through 
space. Even though this postulate is known to be unrealistic in terms of the 
quantised, discontinuous nature of matter, its usefulness lies in the 
simplicity with which macroscopic phenomena may be described. The following 
statement of the first postulate is adopted: 

Postulate I Matter is continuous and distributed uniformly within an 

arbitrary bounded space. 

This statement is more restrictive than those usually offered (S181, 
ZH76) since the uniform and continuous distribution of matter is postulated 
only within a space delineated by boundaries which is defined herein as a 
discrete volume. Thus a discrete volume admits the existence of 
discontinuities at its boundaries. This means that physical phenomena such as 
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shock waves and numerical phenomena such as volume -averaged property 
discontinuities are accommodated within the piecewise continuum model 
postulated. 

Having thus defined discrete volume, the essential requirement is to 
describe the temporal variat of intensive properties within the discrete 
volume from macroscopically oos^ vable conditions. 

For a generalised scalar, ctor or tensor quantity 0 defined by: 


- iKx,t) 


( 2 . 1 ) 


the total temporal derivative of for a cohesive material body is given by: 



* 

• 

ddt 

ipdV = 

(drji/dt) dV + 


v <.m) 

v (m) 


V>(v*n)d A 

( m ) 


( 2 . 2 ) 


Equation (2.2) indicates that the total change of ^ for the entire 
material body is a function of the change of \f) at each fixed point within the 
body plus the transport of at the boundaries of the body. This equation, 
which is known as the 'transport theorem' ( S 18 1 ) places no restrictions on the 
nature of the body other than it be regarded as an autonomous entity within a 
given discrete volume and that it have the characteristics of a 

continuum. In particular, the degree of cohesiveness of the body is 
arbitrary, so that generalisation to a system of particles of arbitrary 
cohesiveness (that is, liquids, solid s or gases) yields: 


ddt 

dV = 

(difr/dt) dV + 

lK v (s>*n)dA 

(2.3) 


V(s> 

V(S, 

A (S) 



Equation (2.3) is known as the 'generalised transport theorem' (S181) and, in 
essence, is the symbolic realisation of the first postulate. 

The transport theorem of equation (2.2) provides the means by which 
macroscopic conservation postulates may first be transformed into their 
microscopic or differential counterparts that apply within the discrete 
volume. Thereafter, the generalised transport theorem permits the 
differential conservation balances to be applied to a system of particles such 
as that comprising the working fluid of a Stirling cycle machine. 

The conservation postulates are expressed strictly in terms of 
macroscopic phenomena. Hence the macroscopic conservation of mass for an 
arbitrary material body is expressed by the following postulate: 
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Postulate II The mass 

time . 

Symbolically, this 


of an autonomous material body is independent of 


may be expressed as: 



(2.4) 


Choosing \j> - p (that is, mass per unit volume or density is the 
transport property) in equation (2.2) produces the differential conservation 
balance : 


dp/dt - - (V* pv) 


(2.5) 


in a Eulerian frame of reference. 

The integral mass balance applicable to a discrete volume is obtained 
from equations (2.5) and (2.3) and is given by: 

dM (s) / dt - p { (v-v (s) ) • -n}dA (2.6) 

J^CS) 

In this generalised or combined Eulerian/Lagrangian form, the rate of 
change of mass of a system of particles is equal to the net advection of mass 
across the boundaries of the system. The advection velocity is the relative 
velocity between the particles and the boundary itself. 

A statement of the macroscopic conservation of momentum for an 
arbitrary material body expresses the third postulate, which is generally 
referred to as Euler's first law (S181) : 

Postulate III The time rate of change of the momentum of an autonomous 

material body relative to the fixed stars is equal to the sum 
of the forces acting on the body. 

Postulate III may be expressed symbolically as: 


ddt 


pv dV — 

* 

sdLA + 

V (W) 

A W) 


V 


pfdV 
( m ) 


(2.7) 
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where s denotes the contact forces per unit area and f denotes the external 
and mutual forces . 

Selecting \j) - pv (the momentum per unit volume) in the transport 
theorem (equation (2.2)), allows the following Eulerian differential momentum 
conservation equation to be derived: 


A 

d(pv)/dt + V*(pw) — V*T - VP + pf 


( 2 . 8 ) 


It should be noted that the sign convention adopted for the extra shear stress 
tensor T is such that T represents the stress acting at any point within a 
material body. Modifying Stokes' hypothesis (Sc79) for gaseous fluids by 
including a 'bulk viscosity' (BS60) allows T to be expressed by: 


T - /i{ Vv + (Vv) T ) + { (A-2^/3)(V.v) }I 


(2.9) 


Together, equations (2.8) and (2.9) represent the Navier- Stokes equations, 
which are usually independently derived on a more intuitive basis (Sc79, 
BS60) . 


The generalised combined Eulerian/Lagrangian form of the integral 
momentum conservation balance is obtained by substituting equation (2.8) into 
equation (2.3) which produces: 


d( [ V ]VM( 5 ))/dt = 


r 


r 


pv{ (v-v (s) ) • -n}cL4 - 

PndLA - 

•M(S) J 

^(S) 1 


A 


+ 

pfdV 



Vis> 



(T* -n)dLA 

<S) 


( 2 . 10 ) 


Thus the rate of change of momentum of a system of particles is equal 
to the net advection of momentum across the boundaries of the system relative 
to the system boundary velocity plus the contact, mutual, and external forces 
acting on the system. 

The fourth postulate is defined by the conservation of energy for a 
material body. In the context of a discrete volume analysis, the first law of 
thermodynamics which is adopted by most authors as their postulational basis 
(Sc79, for example) is not specific enough for a macroscopic material body 
(Go87). Thus the following formulation advocated by Slattery (S181) is 
preferred: 
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Postulate IV 


The time rate of change of the internal and kinetic energy of 
an autonomous material body relative to the fixed stars is 
the sum of the rate at which forces acting on the body do 
work on the body and the rate of energy transmission to the 
body . 

Symbolically, this may be expressed as: 


p(U+v 2 /2)dV = 

(vs)dA + 

A 

p(v*f)dU + 

qdA + 

A 

pEdV 

(2.11) 


A m ) 


A (.m) 

Vim 



where s and f are defined for equation (2.7), q denotes the contact energy 

A 

transmission rate per unit area, and E denotes the external and mutual energy 
transmission rate. The first and second terms on the right hand side 
represent the work done by the corresponding force terms in equation (2.7) 

Choosing the internal plus kinetic energy per unit volume as the 
transport property by setting (V> - p(U+v 2 /2) ) in the transport theorem 
(equation (2.2)), the following Lagrangian differential energy equation may be 
derived (Go87): 


pD(U+v 2 /2)/Dt = p{(v*f)+fb + V* (T*v) - V-(Pv) - V-q 


( 2 . 12 ) 


This equation describes the differential conservation of thermal and 
mechanical energy. It can be simplified by observing that the differential 
conservation of mechanical energy may be determined separately using postulate 
III. Forming the scalar product of the Lagrangian form of equation (2.8) with 
v yields the differential conservation of mechanical energy equation in a 
Lagrangian frame of reference: 


A 

pD(v 2 / 2)/Dt = V-(T-v) - (T:Vv) - (v«VP) + p(v* f) 


(2.13) 


The second term on the right-hand side is a tensor scalar product which 
represents the irreversible conversion of mechanical energy into thermal 
energy, or dissipation. Subtracting equation (2.13) from equation (2.12) and 
expressing the result in Eulerian terms produces: 


A A A 

d(pU)dt + V'(pUv) = pE + (T:Vv) - P(V»v) - V*q 


(2.14) 


This equation defines the differential conservation of thermal energy. It may 
be noted that equations (2.12) and (2.8) contain a redundancy that is absent 
from equations (2.14) and (2.8). Although either set is admissible, and both 
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sets must ultimately yield identical results, the latter set is preferred 
because of the computational simplicity and convenience it affords in 
describing Stirling machine fluid dynamics. 

Choosing the internal energy per unit volume as the transport property 
in the generalized transport theorem (that is, yj> = pU in equation (2.3)) and 
substituting equation (2.14) produces the general Eulerian/Lagrangian integral 
form : 


d( m UM iS) )/dt 


l P E + (T : Vv) + (v-VP)} dV + (q*-n)cL4 

S) J A (S) 


+ 


pH{ (v-v (S) ) - -n)dA 
J A (S) 


? ( V (S)* -n)dA 

. A (S) 


(2.15) 


Thus the rate of change of internal energy of a system of particles is equal 
to the sum of the following components, which in left to right sequence are: 

the rate of mutual and external energy transmission to the system 
the rate of irreversible conversion of mechanical into thermal 
energy within the system boundaries 

the isentropic heat generation rate within the system boundaries 
the net rate of contact energy transmission across the boundaries 
of the system 

the net advection of enthalpy across the boundaries of the system 
relative to the boundary velocity 

the net rate at which mechanical work is done at the boundaries of 
the system. 

In order to implement the integral balances for gaseous fluids, an 
equation of state is required. In keeping with the established practice for 
Stirling machine analysis (Scl871, Wa73) , the ideal gas equation of state is 
used here, namely: 


PV = MRT 


(2.16) 


Nevertheless, there are no intrinsic restrictions placed on the form of 
the equation of state; other equations describing the behavior of real gases, 
such as that of Redlich and Kwong (RK49) , may be used. Owing to their 
relative complexity, however, such equations are not as numerically convenient 
as the ideal gas equation. 

Equations (2.6), (2.10), (2.15), and (2.16) thus provide the analytic 
basis in terms of a discrete volume continuum model for determining the 
working fluid behavior of Stirling cycle machines. 
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2.3 


THE TURBULENCE MODEL 


The integral balances of equations (2.6), (2.10), and (2.15) are 
strictly applicable in the limit as At -*• 0 (Hi75). However, when the balances 
are applied to systems in which At is finite, then the balances are precise 
only for laminar flow conditions. Under turbulent flow conditions, the 
transport properties xj> may experience random fluctuations with periods less 
than At, thus invalidating the instantaneous constancy of the temporal 
gradients implied by the equations as described. The instantaneous value of xj> 
may be represented as the sum of a time-averaged component and a fluctuating 
component : 


xj) — rf> + xf>' 


(2.17) 


Attention here is focussed on obtaining the time-averaged quantities xj> 
directly since the computational effort necessary to obtain V is currently 
beyond the scope of practical Stirling machine simulation. The most general 
approach to obtaining the time -averaged or turbulent transport balances is to 
perform the averaging process on the integral balances directly (S181) . This 
admits fluctuating discrete volume geometries such as those occurring in 
Lagrangian systems. Time -averaging equations (2.6), (2.10), and (2.15) 
results in: 


Mass : 


dM (s) /dt ~ 


Pi (v-v (S) ) • -n)dA 

J^S) 


(2.18) 


Momentum: 


d( IV ]vM (S) )/dt 


h 

i 

( 

pv{ (v-v (S) ) • -n}dA - 

PndA - 

(T* -n)dA 

J^CS) 

A (S) 

A (S) 


p fdV 



(2.19) 


10 



Energy : 


d( [ V ]UW (S ))/dt - {p£ , + (T:Vv) + (v*VP)dU - 

. v t8) 

A 

+ pH{ (v-v (5) ).-n}dA - I P(v (s) *-n)dA (2.20) 

J^CS) J^(S) 

Equations (2.18) through (2.20) are by definition also applicable under 
laminar flow conditions, since from equation (2.17) ^ when “ 0. 

The principal difficulty in solving the time-averaged integral balances 
is the unavailability of the fluctuating component of the transport 
properties. In order to obtain the time -averaged properties directly, the 
analytic approach adopted for dealing with the unknown fluctuating components 
falls into the category of 'Reynolds averaged equations' (Fe83) . In this 
context, since the integral balances derived include both time- and volume- 
averaging, Ferziger makes the following observation: 

' • * • The equations describing the mean field contain averages of 
products of fluctuating velocities and there are fewer equations than 
unknowns - the well-known closure problem. In fact, the set of 
equations can never be closed by further averaging; a closure 
assumption, or what is the same thing, a turbulence model, has to be 
introduced. The closure assumption must represent the unknown higher- 
order average quantities in terms of lower-order quantities that are 
computed explicitly. ' 

The minimum set of assumptions constituting the turbulence model 
adopted here is stated in terms of the six restrictions discussed below. 

These restrictions enable a time-averaged solution of the integral balances to 
be obtained with reasonable computing resources. It must be emphasized, 
however, that the turbulence model adopted is not definitive and is subject to 
amendment by experimental data. 

Restriction I 

The turbulent flow is stationary such that: 

A t turbulence characteristic <K analytic time increment (2.21) 

Restriction II 

The turbulent flow field is spatially homogeneous such that: 
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( 2 . 22 ) 


AV turbulence characteristic <K AV discrete volume 

Rigorously, equations (2.21) and (2.22) are conflicting conditions for 
any turbulent flow because if such a flow field is homogeneous , then it is 
simultaneously a decaying flow field. However, if it is also stationary, then 
the dissipation in the field can only be balanced by a non-homogeneity in 
order to maintain the decaying characteristic. The following rationale 
offered by Hinze (Hi75) is adopted here for proceeding with the stationary, 
homogeneous flow field model: 

'... Fortunately, the rate of decay of the mean properties is rather 
slow with respect to the time scale of the smaller eddies. Therefore, 
the actual state of non-stationarity is considered not to be a serious 
drawback in the experimental study of the smaller scale turbulence. 

For the theoretical study, this makes it possible to apply the concepts 
and theories of stationary random processes.' 

Restriction III 

The stationarity of the turbulent flow field is sufficient for the 
equality of the first and second order time averages, or equation 
(2.21) implies that: 


i> - $ 


(2.23) 


Hence by taking the time average of equation (2.17) it immediately follows 
that: 


- 0 


(2.24) 


Restriction IV 

The ergodic hypothesis is valid for scalar turbulent fields. 

The ergodic hypothesis states that for a stationary and homogeneous 
turbulence : 


mtf* 


V* ” [ensemble]^ 


(2.25) 


where t/>* represents a scalar or component of a vector. Hence from restriction 
IV and equ<_ “ions (2.23) through (2.25) it follows that: 


[VjV’scalar “ ® 

Hence, in particular, for density and temperature: 


(2.26) 
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[V ] P * “ [ V]T 1 


0 


(2.27) 


Equations (2.27) constitute two of the explicit restrictions placed on 
the simulation by the turbulence model, namely, that the temporal fluctuation 
of the volume -averaged density and temperature are zero under the restrictions 
of equations (2.21) and (2.22). Since the ergodic hypothesis is applied to 
scalar properties only and equation (2.26) is limited to volume averages, non- 
zero fluctuations of volume -averaged vector fields and correlations involving 
non-volume- averaged scalar fields are permissible. This requires another 
restriction, namely: 

Restriction V 


The effect on the mean flow resultant from vector turbulent fields may 
be modelled. 

In particular, the effect of the turbulent velocity field fluctuations on the 
time - averaged flow field may be determined using a model such as an empirical 
correlation (for example, that between friction factor and Reynolds number) or 
a two -parameter turbulence model of the k-e variety. 

The last restriction is defined by: 

Restriction VI 


The discrete volume boundaries do not experience temporal fluctuations. 
This may be expressed symbolically as: 


A' 


(S) ~ 


0 


(2.28.1) 


which in turn implies that: 


V* (S) - 0 and 


(2.28.2) 


Equations (2.28) do not place a restriction on the Lagrangian condition 
v ” v (s) as suggested by equation (2.3). Under turbulent conditions, this is 
achieved by setting v (s) - v and admitting a turbulent flux as a function of 
v " v (s> (which is equal to v') across the Lagrangian boundary. In effect 
this converts a turbulent Lagrangian boundary Into a combined 
Eulerian/Lagrangian boundary. In the context of Stirling machine numerical 
analysis, restriction VI enables generality to be maintained without the 
necessity of allowing turbulent discrete volume boundary movement, enabling a 
significant simplification of the subsequent analysis. 
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At this stage , it is convenient to introduce the mass flux g into the 
equations as a means of simplifying the numerical model, g is defined by. 


g - pv 


(2.29) 


Applying equation (2.29) and the six restrictions of the turbulence 
model to the time -averaged conservation balances of equations (2.18) to (2.20) 
yields (in combined Euler ian/Lagrangian form): 

Mass : 


dM (s) /dt 


J ( (g'P v ( sj) • -n)dA 

A (.S) 


Momentum : 


d([tv ]8 ^cs))/dt “ 


l_ g{ (v-v (S) )* -n)dA - _ PndA 

J A (S) J A (S) 

{ (T+T (t) ) • -n)dA + |_ pfdV 

. A ( S) Jy<s) 


(2.30) 


(2.31) 


Energy : 


W (S) )/dt - 


{pE + (T: Vv) + (vVP))dV + 
.^CS) 


_ { (q+q (t) ) • -n)dA 

J A (S) 


_ (Pv (s) *n)dA + 

J^S) 


* 

C P _ (r(g-pv (s) )--n)dA 

J' A (S) 


(2.32) 


Expressing the equation of state (equation (2.16)) in volume -averaged 
terms and applying the turbulence model yields : 


ttV]^ “ [tV]P^[tV]^ 


(2.33) 


The Reynolds stress tensor T (t) in equation (2.31) is given by: 


T (t) - (g v - g v (S) ) - (gv - gv (S) ) 


(2.34) 
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Analogously, the turbulent enthalpy flux vector q (t) in equation (2.32) is 
defined by: 


A __ _ _____ 

q (t) - CpiigT - pv (s) T) - (g T - pv (s) T )} 


(2.35) 


Equations (2.34) and (2.35) cannot be solved since they contain 
additional unknowns for which there are no additional equations (the closure 
problem). Hence a turbulence model is required to solve these equations. 


2.4 THE SPATIAL DISCRETISATION SCHEME 

In his comprehensive review of computational fluid mechanics, Roache 
(Ro82) describes a variety of spatial and temporal discretisation schemes for 
the numerical application of the differential conservation balances which may 
be loosely categorised into 'coincident' and 'staggered' mesh systems. These 
mesh systems may be applied both spatially and temporally, thus yielding a 
multiplicity of schemes involving coincident and staggered mesh systems. 

In a coincident mesh scheme, the three differential transport 
properties- -density , velocity, and temperature- -are evaluated at the same time 
and/or at the same spatial location. However, in a spatially staggered mesh, 
generally the density and temperature are computed at one set of grid points 
while the velocities are computed at an offset grid point mesh. In a 
temporally staggered mesh, the velocity is computed at a half time step offset 
from the density and temperature. In recent years, a de facto consensus has 
emerged that spatially staggered, temporally coincident discretisation schemes 
are convenient and useful for fluid flow modelling (Fe83, Pa80) . Roache 
infers that the first use of a version of the spatially staggered mesh may be 
attributed to Harlow and Fromm (HF64) although its apparent reinvention over 
the intervening two decades is an attestation of its efficacy. 

In the field of Stirling machine analysis, Urieli (Ur77) applied the 
temporally coincident, spatially staggered grid to a simple generic machine 
geometry. Schock, by contrast, used a temporally and spatially coincident 
grid and invoked volumetrically weighted averages for computing flow rates 
between grid points (Sc78) . Both of these discretisation schemes involve the 
application of the differential conservation balances in a one-dimensional 
Eulerian frame of reference. 

Although the spatial discretisation scheme used in the simulation model 
contains elements employed by Urieli and Schock, it has Its origins in the 
mesh structure used in the 'Marker and Cell' (MAC) method of Harlow and Welch 
(HW65) . In particular, the general precepts of a temporally coincident, 
spatially staggered numerical discretisation scheme suitable for the 
application of the differential conservation balances are applied to the time- 
averaged integral balances described In section 2.3. This application admits 
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a three-dimensional space in a combined Eulerian/Lagrangian frame of 
reference . 

As the mass, momentum, and energy balances are fundamentally based upon 
the concept of a discrete volume, the spatial discretisation scheme naturally 
devolves to partitioning a given space into an assemblage of discrete volumes 
with coincident boundaries. All the intensive parameters are thus expressed 
in volume -averaged terms so that, by definition, the value of any intensive 
parameter at a point within a discrete volume is extrinsic to the conservation 
equations. Estimates of the intra-discrete volume intensive parameter 
distribution may be made using 'volume functions' whose complexity is a 
function of the boundary conditions imposed on the discrete volume. 

The essence of the spatial discretisation scheme involves a method of 
constructing the discrete volume grid so that, within any discrete volume, 
scalar intensive parameters are assigned a position and vector intensive 
parameters are assigned a plane, respectively. This introduces the concept of 
a 'volumetric filtration' process which may be defined for a scalar field as 
the assignation of a point to a discrete volume so that volume -averaged scalar 
intensive parameter at that point satisfies the mean value theorem (TM72) for 
the discrete volume. Similarly, for a vector field, volumetric filtration 
assigns a plane to a discrete volume such that the value of a vector intensive 
parameter over the plane (and perpendicular to it) normalised with respect to 
the area of the plane satisfies the mean value theorem for the discrete 
volume . 


The characteristics of the discretisation scheme may be illustrated by 
a two-dimensional Eulerian space using a triangular mesh as shown in figure 
2.1. The mass and energy integral balances are applied to a common discrete 
volume while the momentum balance is applied to discrete volumes which are 
offset from the mass/energy discrete volume and straddle its boundaries. In 
this arrangement, the momentum conservation balance is applied to the net 
momentum, that is, to as many components as the dimensionality of the field 
requires (two in this case). Typically, it is convenient to resolve the 
momentum into components perpendicular and parallel to the plane generated by 
the volumetric filter. 

However, although generalised, this scheme is not necessarily the most 
efficient numerically. Consider, for example, the Cartesian mesh shown in 
figure 2.2 (the same observations apply to any regular mesh such as 
cylindrical or spherical). In this case, the net momentum balance may be 
split into its vector components so that each component balance is applied to 
a unique offset discrete volume. As indicated by the double cross-hatched 
area in figure 2.2, the net momentum can be calculated by the vector addition 
of the momentum components determined individually for that volume. Thus this 
particular case still yields the net momentum over the entire mass/energy 
discrete volume as does the general scheme but requires half the computational 
effort in two dimensions. 
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Figure 2.1 Spatial discretisation scheme characteristics. 


X 2 




MASS/ENERGY DISCRETE 
VOLUME 



MOMENTUM DISCRETE 
VOLUME IN THE 
x DIRECTION 



MOMENTUM DISCRETE 
VOLUME IN THE 
x DIRECTION 


x 


1 


Figure 2.2 Cartesian mesh spatial discretisation. 







Generalising these concepts to a three-dimensional space with an 
arbitrary coordinate system, the characteristics of the staggered discrete 
volume grid may be expressed symbolically by denoting the mass/energy discrete 

volume as V {S) and any momentum discrete volume as V n(S) . Then considering any 
adjacent mass/energy discrete volumes 17 {S) 2 and V (S )i+i separated by a boundary 
A rus)i + 1 ( w ^ ere represents an arbitrary sequence) , the following attributes 
may be defined for the discrete volume grid. 

Attribute I 

A momentum discrete volume straddles every mass/energy discrete volume 
boundary A n(S)i+1 and is bounded by a convex surface containing the 
perimeter of the boundary ^n(s)i+i an< ^ the volume centroids of the 
adjacent mass/energy discrete volumes. Symbolically this becomes: 

[ a < v n<s>i+l> n < J <W> u 3 < U <s )i+i) ) 1 c u ^(^<sji+i)l (2.36) 

Attribute II 

The volume -averaged intensive parameters corresponding to a mass/energy 
discrete volume are located at its volume centroid. Thus [V ]V> is 
located at tV ]X such that: 

[V] x “ ( xdV) / V (2.37) 

. V 

Attribute III 

The volume -averaged intensive parameters corresponding to a momentum 
discrete volume are located on the boundary d n(S) ^ +1 . Therefore iv n ] V* is 
determined by: 

f i/>dV - if>dV + ifrdV (2.38) 

^n(S) 

where a,p < 1 satisfy equation (2.36). 
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Attribute IV 


The net momentum of any mass/energy discrete volume is uniquely defined 
everywhere in the volume by the momenta of its boundary momentum 
discrete volumes, or: 


Sn^ncs) (2.39.1) 

r ms> 

where : 


S a <^n(s)) “ *<?(«>) (2.39.2) 


. 1 


gdV 


CS) 


Applying these attributes to the particular Cartesian mesh of figure 
2.2 simplifies equation (2.36) to: 


* _1 ( *OW> n HV ms)i+l ) ) - V isU /2 (2.40.1) 

^(s)i + ^<s)i+i “ v ms ) (2.40.2) 


which satisfies equations (2.37) and (2.39). 

Equation (2.37) implies that in the case of the one-dimensional 
Cartesian coordinate system which is most often used for Stirling machine 
analysis, a mass/energy discrete volume centroid is defined by a plane 
separating two adjacent discrete momentum volumes. In a two-dimensional, 
Cartesian system the centroid becomes a line perpendicular to the mesh surface 
while only in three dimensions does the centroid become a point. 

The attributes of the spatial discretisation scheme enable all the non- 
volume-averaged boundary terms in the time -averaged integral balances of 
equations (2.30) through (2.32) to be replaced with volume -averaged terms, so 
permitting the balances to be numerically discretised directly. The final 
forms of the integral balances used in the simulation model are given in 
combined Euler ian/Lagrangian form by: 

Mass : 


d M is) /dt 


_ {([tv n ]g - [tv n ]P v n(S) ) . -n)dA 

J^ms) 


(2.41) 
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where the subscript n denotes that the relevant parameters are associated with 
the momentum discrete volume. 

Momentum : 


d([ t v n ]g v rus))/ dt 


_ [tv]g{ ( [tv] v ' v (S)> * -n)dA 

J A (.S) 


I [tVJ PnciA 

J^CS) 


- _ ( m T (t) *-n)dA - W n(S) g 

J A (S) 


(2.42.1) 



( [ V [tV] v + ( V [tV3 v ) T ) 


+ l ((tv]A~2 [tV j/i/3) (V* [t v]V) )I] • -n 


dA 


where T in equation (2.31) is replaced by the volume and time average of 
equation (2.9) and the Reynolds stress tensor is given by the correlation: 


m T(t> " -m8' m v ' 


(2.42.2) 


Energy: 


C v d( ltV] T M (S) )/dt - V (s) { [tv] pf + [tv] (T:Vv) + 


+ ([tv] 


— [tv n ] K (^[tv n i^* ~ n ) 

J ^ms) 



{( ( v n ]q <t) * 

( S ) 


-n)dA 


(2.43.1) 


A 

+ Cp 


_ [tv n ) r Httv n ]g- [tv n ]P v n(s)^ 
•Mn(S) 


-n)cLA 
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where defines the turbulent dissipation tensor which is modelled and 

the turbulent flux vector takes the form: 


A — — - — — 

[V n ]qCt) “ C p(lV n )P [V n jV' [V n ]^ - . ]P [V n ]V' [ t V n ]^) 


(2.43.2) 


2.5 THE NUMERICAL SOLUTION ALGOR THM 

An implicit (or advanced time) numerical algorithm is used to 
temporally integrate the total time-differential conservation balances of 
equations (2.41) through (2.43). The algorithm has its origins in the 
Implicit Continuous -fluid Eulerian (ICE) technique of Harlow and Amsden (HA71) 
and the Semi - Implicit Method for Pressure -Linked Equations, Revised (SIMPLER) 
of Patankar (Pa80) . 

Central to the algorithm is the notion (HA71) that a change in the 
pressure field is a function of the information propagation rate, or: 


A P = (dP/d p)jAp 


(2.44) 


where ( dP/dp) T is the square of the isothermal speed of sound for a fluid with 

constant specific heats. Substituting the equation of state (equation (2.33)) 
and discretising yields: 

V / R [t v]T - M s ) / A t « d M (S )/dt (2.45) 

after multiplying through by V / At . Equating with equation (2.41) produces: 


rtv]^ V/R ltV] T At “ M s /At 4- 


<<ttv n ]g 
J A IUS) 


[t v n ]P v n(S)^*' n ^ dA 


(2.46) 


Equation (2.42.1) may be discretised and rearranged to produce: 


lt V g 


([tv n ]g v n<.s)) a / v rus) 


At{ 


_ [tvjPndA + f( (t v ]g. 

Ks) 


[tvj^ . 


imP ) ) / v ms) 


21 



Substituting into equation (2.46) for [ tv ^]g and rearranging: 


ttvi* V / R [tvl r At + (A t/V ms) ) 


rus) 


(i 


_ itviPndA • -n 


dA 


Hs) 


M a /At + 

. ^n(S) 


Equation (2.47) yields an advanced time or implicit pressure field 
which is obtained by linking the pressure terms in the conservation of mass 
and momentum balances, hence the 'pressure-linked' terminology. However, 
equation (2.47) can only be solved given the advanced time temperature and 
density fields which in turn depend primarily on the advanced time mass flux 
field. This dilemma is resolved by using an iterative solution scheme based 

on an estimated or guessed mass flux field ftv n i8* as used i n the SIMPLER 
algorithm. 

The set of steps comprising the algorithm may be described broadly as 
follows (see Go87 for specifics): 

1. Guess the mass flux field [ t v n ]8*- 

2. Explicitly compute the discrete volume masses from equation (2.41) 
and then infer ttV] p from known values of I7 (s) . 

3. Implicitly compute the temperature field using equations (2.43). 

4. Implicitly compute the pressure field from equation (2.47). 

5. Implicitly compute the mass flux field from equations (2.42) using 
the pressure field determined in step 4. 

6. Compare [tv jg (computed) with [tv n )8 (guessed) and return to step 2 

with f( t tv ]8.[tv n ]8*) -*• ttv n ]8* i f the mass flux field is 
insufficiently converged. 

As with all implicit schemes, this algorithm requires the repeated 
numerical inversion of matrices. Hence the cost-effectiveness of the 
algorithm for transient compressible fluid flow simulation is limited by the 


[l([tv n ]8 v ms ))“ + At f([tv n j8. [tvi r » ttvjP) ) / v ms) 


tty 


,\P v ncs>] * ~ n 


dA 


(2.47) 
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size of the matrices generated. This arises because a break-even point is 
eventually reached when the cost of matrix inversion becomes equal to the cost 
of using an explicit algorithm with time steps small enough to satisfy the 
Courant criterion (CF67) . However, for all the problems described in this 
report, the algorithm is constrained to operate within the cost-effectiveness 
break-even point by limiting the spatial discretisation (particularly in the 
two-dimensional problems) . 

It should be noted that equations (2.45) and (2.47) are discretised 
using a simple first order temporal difference. In general, this is not 
necessary and other higher order temporal discretisation schemes may be used 
(see section 4.7). 


2.6 INFORMATION PROPAGATION MODELS 

The viability of an implicit analysis for compressible fluid flows 
depends on the extent to which the analysis properly accounts for information 
propagation both on a cyclic equilibrium as well as on a transient basis. In 
an explicit analysis, accurate information propagation modelling requires that 
the Courant criterion (CF67) is met at each discrete volume. This imposes 
limitations on the integration time increment for a given spatial 
discretisation. The issue then is to find a similar criterion for selecting 
the time increment in an implicit analysis. Consider equation (2.47) in the 
following reduced form: 


i+m 

E K j ltV}Pj 
j—i-m 


(2.48) 


where i denotes the individual mass/energy discrete volumes of which the flow 
area is comprised, m is dependent on the dimensionality of the problem, and K • 


and depend on At. On a transient basis for arbitrary At, equation (2.48) 
can be applied to a series of pressure domains, one for each discrete 
mass/energy volume as illustrated (for a particular two-dimensional Eulerian 
field) in figure 2.3. 


J 


Each pressure domain has an extent ?f(v±v a I)^At determined by the 
information propagation characteristics where (v a ) ± is the sonic velocity 
within each discrete mass/energy volume comprising a particular pressure 
domain. Hence by this process of partitioning, a 'pressure domain splitting' 
(PDS) algorithm may be structured to model information propagation phenomena. 
This simplified explanation ignores the complexities arising from defining the 
pressure domain boundaries under supersonic or sonic flow conditions (that is, 
when | v | ^ > (v a ) . However, it may be mentioned that under these conditions 
the PDS algorithm essentially devolves to a standard approach such as the 
9 region- to-region' method (Jo69) . Unfortunately, the PDS approach is 
computationally quite expensive and hence may not be practical for Stirling 
machine analysis. 
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Figure 2.3 Structure of the pressure domain splitting algorithm. 


However, the principal focus of Stirling machine simulation is the 
portrayal of cyclic equilibrium conditions rather than the details of a 
specific transient cycle. This suggests a significant simplification by 
applying equation (2.43) to the entire flow field (or treating the flow field 
as a unitary pressure domain) . Two approaches to obtaining the cyclic 
equilibrium solution under these conditions may be hypothesized: 

1. Infinite information propagation. 

This hypothesis states that, at cyclic equilibrium, sufficient 
time has passed such that every point in the flow field has received 
information from every other point in the flow field for all instants 
over the cyclic period. The hypothesis may be implemented by 
arbitrarily selecting an integration time increment that is much less 
than the smallest information propagation time characteristic of a 
particular machine. The time characteristic may be defined as the 
interval required for a pressure wave to exactly traverse the unitary 
pressure domain once. Henceforward, the infinite information 
propagation hypothesis is distinguished by referring to its 
implementation as the 'equilibrium algorithm . 

2. Characteristically determined integration time increment. 

Here the integration time increment is treated as a dependent 

variable which is instantaneously equal to the machine time 
characteristic. This yields the spatially limiting case of the PDS 
algorithm and as such approximates the full information propagation 
transient solution. The implementation of the characteristically 
determined integration time increment hypothesis is termed the 'unitary 
pressure domain' (UPD) algorithm. 
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The characteristic number N c ^ is a convenient way of describing the 
information propagation characteristics of a Stirling machine. N ch is defined 
as the number of complete pressure traverses between the expansion and 
compression space piston faces occurring over a cycle or: 


N ch “ 1 / f ?[ 2 i / min { \ v i + ( v ah I > | v i " ( v a) i | ) ] (2.49) 

Hence the larger N c ^ , the more accurate the equilibrium algorithm is likely to 
be in predicting cyclic equilibrium performance. 

2.7 BOUNDARY ADVECTION 

In the light of the attributes of the spatial discretisation scheme 
and, in particular, the area normalisation requirement for the volume - averaged 
momentum, the interpretation of the momentum advection term in equation 
(2.42.1) is not self-evident. The usual approaches to boundary advection do 
not translate directly into the discrete volume integral framework where 
boundary discontinuities are admissible although not necessary. Several 
approaches to dealing with boundary advection are possible. These include a 
simple one-dimensional equilibrium approach, a multi-dimensional volume 
function analysis and flux corrected transport methods such as those of Book, 
Boris and Hain (BB75) and MacCormack (Mc82) . In the simulation model 
described in this report, the simple one - dimensional equilibrium approach has 
been adopted as a baseline and modified by empirical models where necessary. 

Physical insight to the form of the boundary advection of momentum may 
be obtained by solving equation (2.42.1) under equilibrium conditions ignoring 
mutual and external forces and turbulent momentum fluxes, that is: 


_ [tv ]g([tv ]V*-n)dA - |_ ( ttv } T*-n)dL4 = 0 

ks, n 


(2.50) 


Considering the pair of adjacent discrete momentum volumes shown in figure 2.4 
in one dimension, equation (2.50) has the general solution (dropping the 
averaging notation for clarity) : 


(g A ) ~ KgA) nR -(gA) nL ){exv(N pe x/l) - 1) / ( exp(N pe ) - 1 ) + ( gA) nL (2.51.1) 


25 



I 



Figure 2.4 Adjacent one -dimensional momentum discrete volumes. 


where Np e is the Peclet number at the centroid of the discrete mass/energy 
volume of length I given by: 


N Pe - 3pv2/4/i 


(2.51.2) 


Following Patankar (Pa80) and denoting the net boundary momentum flow 
(advection plus diffusion) across boundary A in figure 2.4 as G , equations 
(2.51) produce: 


if N Pe - 0 then: 

G - (4 M /3 pD((gA) nL -(gA) nR ) 

if N Pe i 0 then: 

G - v[(gA) nL +{(gA) nL -(gA) nR )(exp(N Pe )-l)] 


These equations provide a physically meaningful methodology for 
determining the boundary advected momentum flux. Consider a plot of G as a 
function of Np e as shown in figure 2.5. Since N Pe expresses the ratio between 
the advection and diffusion of momentum across a discrete momentum volume 
boundary, figure 2.5 shows that even at very low Np e the advection term 
dominates. When Np e - 0 there is no advection while in the intermediate range 
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Figure 2.5 Net momentum flux dependence on Peclet number. 


the momentum flow is partly diffusive and partly advective. Furthermore, in 
the limit: 


As N Pe -* °°> G v (S A ) n L (2.53.1) 

As N Pe G -> v (.S A )nR (2.53.2) 


Thus figure 2.5 provides a model for determining the boundary advection 
of momentum. In keeping with the transient nature of Stirling machine fluid 
flow, the diffusion is kept as a discrete term in equation (2.42.1) and is not 
lumped together with the advection term. Hence for a discrete momentum volume 
boundary contained within a discrete mass/energy volume i, the advected mass 
flux is defined in terms of the relative velocity v^ perpendicular to the 
boundary where : 


v * - [tv]V - v (S) (2.54.1) 

Hence : 
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If > 0 then: 


ttv]gO**-n)dA 
. A (S) 


i 


( v *[tv n )g A n(S))i 


(2.54.2) 


If v^ < 0 then: 



ttv]g( v ** -n)dA 
(S) 


i 


v i([tv ]8 ^n(s)^I+ 1 


(2.54.3) 


where i and 1+1 denote adjacent discrete momentum volumes. 


Equations (2.54) intrinsically perform the area normalisation required 
for equation (2.42.1) as mandated by the attributes of the discretisation 
scheme. This occurs because the transport term is gA (or the mass flow rate) 
which allows [tV ,g to be normalised by the appropriate flow area. 

It may be noted that equations (2.54) represent a convoluted integral 
version of the 'second upwind differencing' method proposed by Gentry , Martin, 
and Daly (GM66) . An analysis of this method shows that while it is 
transportive (as with the classical or first upwind difference), it is also 
second order accurate for the advection field (Ro82) . 

The boundary advection of enthalpy in equation (2.43.1) is determined 
analogously to that described for the momentum equation with a similar end 
result. Thus for a discrete mass/energy volume boundary contained within a 
discrete momentum volume i, the advected enthalpy flux is defined in terms of 
the relative boundary mass flux g ± perpendicular to the boundary where. 


g* “ ([tv n )g ‘ ( tv n jP v n<s>) 


(2.55.1) 


Hence : 

If g L > 0 then: 



[t v n] r(g*-n)dA 

') 


i 


( S * [tv] A n(s)^ i [tvi r i-i 


(2.55.2) 


If < 0 then: 
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(2.55.3) 


_ ttv nJ r(g*-n)dA - (g* [tv ,r A n(s) ) 1 
A ms ) i 


where i-1 and i denote adjacent discrete mass/energy volumes. 

By nature of its simplicity, the approach adopted for modelling 
transient advection does not resolve completely the 'false diffusion' problem 
(Ro82) engendered by any upwind scheme. However, in terms of Stirling machine 
simulation, this problem only becomes predominant for the enthalpy advection 
in the regenerator. Thus in this case, an empirical model is used to overcome 
the false diffusion deficiency of the advection scheme adopted (section 3.5). 
More advanced advection modelling schemes based on volume functions have been 
developed and indeed do yield better results, but in terms of Stirling machine 
simulation, at least in one dimension, it is doubtful whether the improvement 
in accuracy (which is small when the empirical model is correctly adjusted) is 
worth the additional cost of computation. 


2.8 CLOSURE 


The foregoing discussion has described the manner in which the 
fundamental postulates are transformed via a turbulence model and a 
discretisation scheme into equations suitable for direct implementation in a 
computer code. The derivation of the final equations places no restrictions 
on the dimensionality of the implementation- - the same equations are applicable 
to one-, two- and three-dimensional problems. The information propagation 
model and implicit numerical algorithm developed allow the equations to be 
integrated temporally so that the influence of information propagation may be 
bounded. This permits the impact of information propagation effects on the 
performance of Stirling machines to be at least qualitatively investigated. 
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CHAPTER 3 


SIMULATION OF TH E SPACE — P 0 tf E R 

DEMONSTRATOR ENGI N.E 


3.1 


INTRODUCTION 


The primary objective in simulating the Space Power Demonstrator Engine 
(SPDE) was to provide a test bed for validating the turbulence model developed 
from the Mechanical Engineering Test Rig (METR) experimental results The 
SPDE was selected by NASA as the target engine because of its topical 
relevance to their Stirling engine program and also because it is a more 
complex device in terms of fluid dynamics than previous generations of free- 
piston Stirling engines (such as the Sunpower RE1000) . This provides a more 
severe test of the turbulence model and hence increases its potent a 
application envelope. 

The simulation approach adopted has been to model the SPDE as a one 
dimensional system with the option of replacing the one -dimensional heater 
module with a two-dimensional module. The heater was chosen as the two- 
dimensional replacement module since it was targeted for flui ynam c 
similarity with the METR. The purpose of performing a two-dimensional 
simulation of the heater is to determine on a systems basis the effect ot 
using one-dimensional friction factor and heat transfer coefficient 
correlations in a Stirling heat exchanger in comparison with a model that does 
not require such first order empiricism. Nevertheless, the requirement fo 
second order empiricism, in particular the use of turbulence models, remains. 
The development of such a model is described in chapter 5. 

In view of the long lead times expected (and later realised) in 
producing validation-quality METR experimental data upon which an oscillating 
flow turbulence model could be based, prudence dictated. that the SPDE 
simulation codes be developed first using standard friction factor and heat 
transfer coefficient correlations for the one-dimensional components. 
Furthermore, the initial two-dimensional heater module simulations would hot 
include any turbulence modelling, that is, the flow is assumed to be laminar. 
This provides the necessary baseline against which the effect of a turbulence 
model (as well as any improved one -dimensional correlations based thereon) may 

be compared. 

A comparison of one- and two-dimensional heater module simulations also 
provides r basis against which the efficacy of two-dimensional simulation of 
Stirling machines can be judged, that is, whether the improved simulation 
accuracy is worth the increased cost of computation. Such a judgement is no 
Xlous because It can be hazardous to Isolate the effects of oscillating flow 
in Stirling machine heat exchangers on the overall machine performance. This 
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arises since the cyclic energy balance is the product of the synergistic 
interaction of several factors including: 

porous flow in the regenerator, 

compressibility effects, particularly the influence of information 
propagation, and 

multi -dimensional flow fields in the working space cylinders and 
heat exchanger plena. 

Hence it is possible that in a systems context, the advantages of single 
component two-dimensional simulation are nullified by the transport boundary 
conditions imposed by the rest of the closed system. This is an important 
consideration in attempting to apply any correlations experimentally derived 
from the METR directly to Stirling hardware because of the boundary condition 
dissimilarities . 

Thus, the intention of the one- and two-dimensional simulations of the 
SPDE discussed hereafter is to provide an analytic foundation based on actual 
hardware performance for later application of the METR experimental data. The 
validity of this foundation is assessed by a comparison of the baseline one- 
and two-dimensional simulation predictions against the cyclic performance data 
produced by the SPDE test program. Unfortunately, for the reasons noted in 
chapter 1, it did not prove possible to include the turbulence model developed 
in chapter 5 into the SPDE two-dimensional simulation, nor was it possible to 
go one step further and develop one -dimensional correlations from either the 
simulated or analytic turbulence data. 


3.2 SIMULATION HARDWARE 

A schematic of the SPDE is shown in figure 3.1. The engine consists of 
a pair of back-to-back, beta-configuration, free-piston Stirling engines which 
share a common expansion space. Work is extracted via linear alternators 
attached to the pistons. As the component free-piston engines are 
symmetrical, only half the SPDE need be simulated when using a one -dimensional 
system model even with the inclusion of a two-dimensional heater module. The 
SPDE has a design oscillating frequency of 105 Hz and a pressurization of 150 
bars . 


The expansion and compression spaces of the SPDE are quite complex in 
comparison with those of typical Stirling machines such as the Sunpower RE- 
1000 (Sc83) and GM-GPU3 (Th79) engines. Part of this complexity stems from a 
restriction of the maximum piston and displacer amplitudes to about 12.7 mm. 
The resulting bore- to- stroke ratios predicate that the displacer and piston 
function more as flat-plate oscillators than as conventional pistons. For 
this reason, the flow in the expansion space / heater plenum (figure 3.1) is 
strongly two-dimensional with a significant radial gradient of axial velocity 
at the heater entrance. The compression space also exhibits notable two- 
dimensional momentum boundary conditions in addition to an unusual geometry, 
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Figure 3.1 Schematic of the Space Power Demonstrator Engine. 



namely, two variable volume components joined by a conical annular passage. 
Therefore, in the light of these complications and the importance of correctly 
accounting for two-dimensional momentum boundary conditions (Go87), a strictly 
one-dimensional system model is not a satisfactory description of the SPDE. 

From a fluid dynamics perspective, the balance of the SPDE is of fairly 
conventional design. A tubular heater and cooler are separated by a square 
mesh, woven screen regenerator. The displacer is supported via a gas spring 
while the piston oscillates against a bounce space. Gas bearings are used 
throughout and a bearing gas supply system is built into the engine. As these 
details are not essential for implementing the simulation (which relies on 
predefined displacer and piston harmonic motions), the reader is referred to 
reference Br87 for further information. 


3.3 THE SPDE SYSTEM MODEL 

The overall system model of the simulated SPDE is shown in figure 3.2 
while a listing of the specific geometry used is given in table 3.1. The 
regenerator and cooler are treated as one dimensional while the heater may be 
one or two dimensional. The expansion space is split into a cylinder cavity 
and a heater plenum. The compression space is divided into upper and lower 
variable volume components joined by a conical connecting passage which is 
also given a one -dimensional discretisation. 



exp. space / heater 
plenum 


Figure 3.2 SPDE system simulation model. 
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Table 3.1 SPDE simulation geometry. 


Expansion Space: 


Length 

57.24 mm 

Diameter 

114.4016 mm 

Wetted area 

12.86 cm 2 

Volume at disdacer datum 

1 7 . 1 cm 3 

Heater Plenum: 


Axial length 

It 16 mm 

Axial flow area 

289.44 cm 2 

Axial wetted area 

71.344 cm 2 

Volume 

298.92 cm 3 

Radial length 

54.515 mm 

Radial flow area 

53.944 cm 2 

Radial wetted area 

268.77 cm 2 

Radial entrance area 

36 . 544 cm 2 

Volume within radial momentum 

123.32 cm 3 



Heater: 


No . of tubes 

1632 

Length 

90.17 mm 


1,27 mm 

Regenerator : 


Casing inner diameter 

139.6746 mm 

Casing outer diameter 

223 . 5708 mm 

Length 

25 . 4 mm 

No. of layers in gauze stack 

350 

Gauze mesh 

200 ins” 1 

Gauze wire diameter 

0.04064 mm 

Volumetric porosity 

0.7105 

Gauze density 

7833.03 kg/m 3 

Gauze specific heat capacity 

502.42 J/ke.K 

Cooler : 


No. of tubes 

1584 

Length 

95.25 mm 


1.524 mm 

Upper Compression Space: 


Length 

3.81 mm 

Hydraulic diameter 

21.656 mm 

Midpoint flow area 

238.3 cm 2 

Upstream boundary flow area 

115.95 cm 2 

Wetted area 

41.462 cm 2 

Displacer area 

46.259 cm 2 

Volume at displacer datum 

89.791 cm 3 

Axial momentum volume 

89.791 cm 3 

Flange Passage: 


Length 

24.13 mm 

Hydraulic diameter 

15.457 mm 

Midpoint flow area 

33.042 cm 2 

Upstream boundary flow area 

33.042 cm 2 

Wetted area 

293.71 cm 2 


84 . 102 cm 3 
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Table 3.1 (continued) SPDE simulation geometry. 


Linking Annulus : 



Length 

10.846 

mm 

Hydraulic diameter 

14.566 

mm 

Midpoint flow area 

40.66 cm 2 

Upstream boundary flow area 

32.086 

cm 2 

Wetted area 

121.11 

cm 2 

Volume 

44.101 

cm 3 

Joining Ring Passage: 



Length 

61.488 

mm 

Hydraulic diameter 

11.683 

mm 

Midpoint flow area 

32.086 

cm 2 

Upstream boundary flow area 

32.086 

cm 2 

Wetted area 

656.73 

cm 2 

Volume 

191.81 

cm 3 

Peripheral Compression Space: 



Midpoint flow diameter 

266.31 

mm 

Midpoint flow area at piston datum 

35.156 

cm 2 

Upstream flow diameter 

266.31 

mm 

Upstream flow area at piston datum 

35.156 

cm 2 

Wetted surface diameter 

455.72 

mm 

Wetted surface area at piston datum 

355.58 

cm 2 

Piston area 

108.83 

cm 2 

Volume at piston datum 

169.73 

cm 3 

Central Compression Space: 



Length 

42.385 

mm 

Wetted area 

144.98 

cm 2 

Displacer area 

55.851 

cm 2 

Piston area 

56.438 

cm 2 

Volume at piston and displacer datums 

242.91 

cm 3 


Initially, attempts were made at using a strictly one -dimensional 
description of the upper and lower compression spaces as well as the expansion 
space. However, as expected (Go87) , this approach was not successful in 
enabling the experimental cyclic performance data to be matched. Hence a 
"pseudo- two-dimensional” method of including the actual two-dimensional 
momentum boundary conditions was developed. This method is illustrated in 
figure 3.3. Figure 3.3.1 shows an expanded view of the expansion space 
cylinder cavity and heater plenum. The mass fluxes at planes A and B are 
evaluated using sequential one -dimensional momentum control volumes. However, 
these orthogonal mass flux vectors are advectively and diffusively decoupled 
from each other in compliance with the two-dimensional topology. Similarly, 
the one -dimensional boundary condition within the cylinder for the mass flux 
computed at plane A is g r - 0. Similar principles hold in figure 3.3.2 for 
the upper compression space in which radial and axial velocity and mass flux 
components are maintained in their proper vectorial relationship. The pseudo- 
two- dimensional discretisation of the lower compression space shown in figure 
3.3.3 is achieved by dividing the space into peripheral and central zones. A 
one -dimensional mass flux is thus computed within the lower compression space 
based on physically appropriate boundary conditions. 
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Figure 3.3 


Figure 3.3 


Figure 3.3 


L 


1 Expansion space / heater plenum pseudo- two-dimensional 

discretisation. 



I* 


x 


.2 Upper compression space pseudo- two-dimensional 

discretisation. 



In the two-dimensional heater module, the heater is described as a two- 
dimensional parallel aggregation of the 1632 separate tubes of which it is 
composed. Hence the two-dimensional flow parameters calculated for any 
control volume in the single aggregated tube are assumed to prevail in all the 
tubes. This assumption is rot really satisfactory, in view, for example, of 
the radial gradient of axial locity that exists in the heater plenum for 
positive gas flows (expansion tc compression space) . This implies that the 
flow in the inner rows of tubes i likely to be rather different from those in 
the outer rows of tubes, particul^ ly in terms of turbulence triggering 
effects. However, the better altei lative of specifying several parallel two- 
dimensional heater flow paths was n^ t a pragmatic alternative owing to the 
preliminary nature of this investigation combined with the large computation 
costs involved. Such an approach deserves to be tested in the future. 

The system is discretised spatially by assigning 12, 7, and 3 control 
volumes to the regenerator, cooler, and conical connecting passage, 
respectively. The heater is modelled using 7 axial and 6 radial control 
volumes while all the remaining components are represented by single control 
volumes . 


3.4 SIMULATION CONSIDERATIONS 

The standard set of equations described in section 2.4 is used. In two 
dimensions, the equations are mapped onto a cylindrical coordinate system in 
the heater module. Hence, all the control volume parameters in the two- 
dimensional heater module naturally reflect its inherent radial symmetry as 
shown in figure 3.4. 

Every mass/energy control volume is associated with two axial (x) and 
two radial (r) momentum control volumes in accordance with the discretisation 
attributes listed in section 2.4. A linear radial spatial discretisation 
(constant Ar) was selected because, computationally, it is a more severe test 
of the simulation than other physically more attractive schemes (such as a 
logarithmic radial discretisation with the radius decreasing towards the tube 
wall). This arises because not biasing the radii to discretise more 
accurately the steeper boundary layer velocity gradients stresses the volume - 
averaging procedure inherent in the integral equations more severely, thus 
increasing the likelihood of any errors in the simulation becoming manifest. 
This is particularly true in terms of testing the generality and viability of 
a turbulence model in an integral framework. 
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The radial momentum flux boundary values are defined as follows: 


Tube wall: (g nr ) i Wn - 0 

Axis of symmetry: i 1 “ ® 

LHS boundary: 

if .51 (g nx A nx ) 1 j + (Snx^nx^iJ-i* - 0: 

(8nr)iJ “ 0 

if • 5{ (g nx ^xix^i ,j + < 

(3.1) 

i,j “ 

RHS boundary: 

if . 5{ (gnx^nx^N+i* J + * J " ^ ) < 

(Snr^tf+iJ " 0 

if . 5 { (gnxAnx^N+i > J + ^nx^nx^N+i» — 0: 

^nr^//+i,J " (Snr^N+i,J ^ 


The interface between the one -dimensional and two-dimensional meshes is 
shown in figure 3.5 in terms of an aggregate two-dimensional control volume on 
the expansion side of the heater (the regenerator side interface is a mirror 
image of the one shown) . 



momentum interface 
region 


Figure 3.5 Interface between the one- and two-dimensional meshes. 
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The mass and energy transport interfaces between the one-dimensional control 
volume I and the two-dimensional control volumes X+l,j (where 1+1 «> (i = U) 
are accomplished naturally in terms of the two-dimensional advective fluxes 
(g nx )i j at plane C and the diffusive fluxes defined by the properties at 
planes B and D and the distance between them perpendicular to plane C. This 
arises from the discretisation attributes (section 2.4) which locate the 
volume -averaged intensive properties of control volume X at plane B (as 
opposed to a point, which would only occur, by definition, in a three- 
dimensional discretisation) . The radial momentum interface advective fluxes 
are defined by equations (3.1) and by noting that (9 r nr )j j E 0 in the one- 
dimensional control volume X. Diffusive radial momentum interface fluxes are 
ignored because of their smallness (compared with the advective fluxes) and 
because of the arbitrariness of the definition (g nr )j j - 0 (which is, of 
course, untrue in reality). 

The x momentum interface is more problematic because the definition of 
the advective velocities (v x ) x j in the one-dimensional control volume X is 
somewhat arbitrary. The diffusive flux is computed naturally from the 
gradient between {g nx )i,j at plane C and (g nx ) r at plane A. After numerical 
experimentation, the approach finally adopted for the advective flux interface 
is defined via (v x )j j for the general case ((A nx )j (Ax^X * ( A nx^ I+i^ ^y. 

(v x )j j = (A nx )i,j {(A nx )j i9 nx )x /( A nx)j+i + ( A x^I,j (3. 2) 

This formulation is mass conservative while allowing (v x )jj to vary radially. 
The other basic formulation is also mass conservative but assigns a single 
value to all (v ) r . This approach did not fare as well as the approach 
adopted in tracking the boundary layer growth . 

Standard Kays and London (KL64) friction factor and heat transfer 
correlations are used in all the one -dimensional control volumes (including 
those in the one - dimensional heater module) while all the turbulent terms in 
the two-dimensional equations are zeroed, that is: 


[V] T (t) = [v n) q (t) - - 0 

3.5 RESULTS 


(3.3) 


The two experimental test points used for the validation exercise are 
defined by the parameters listed in table 3.2. Test 46 represents a lower 
power point while test 42 approaches the maximum power output of the SPDE. As 
the test configuration consists of two back-to-back engines, the wall 
temperatures as well as the piston and displacer amplitudes are taken as the 
mean of the left- and right-hand engine parameters. 
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Table 3.2 Simulation input parameters. 


SPDE Test 42 
(11/9/86) 


SPDE Test 46 
(14/3/86) 


Working fluid 

Helium 

Helium 

Frequency (Hz) 

99.385 

99.569 

Charge pressure (bar) 

149.67 

150.29 

Heater wall temperature® (K) 

574.555 

677.89 

Cooler wall temperature® (K) 

308.052 

345.11 

Displacer amplitude® (mm) 

6.927 

7.8016 

Piston amplitude® (mm) 

6.9255 

9.1276 


Notes : 


Taken to be the mean of the left- and right-hand engine test parameters. 


Three simulations were carried out for each test point, UPD and 
equilibrium algorithm simulations for the one -dimensional heater module and an 
equilibrium algorithm simulation for the two-dimensional heater module. The 
non-dimensional parameters and correction factors associated with the 
simulation runs are listed in table 3.3. All the friction factor and heat 
trans fer coefficient multipliers operating on the Kays and London correlations 
are at their baseline values of unity with the exception of those for the 
regenerator matrix friction factor. These latter multipliers reduce the 
nominal steady- state friction factor in the regenerator by 35% and 45% for 
tests 46 and 42, respectively, in order to match the experimental data 
according to the validation protocol developed in Go87. This is thought to be 
a consequence of radially non-uniform mass fluxes in the regenerator (owing to 
entrance effects) and the high frequency of the flow oscillation. 

In view of the large Reynolds numbers resultant from the high frequency 
and pressurisation of the SPDE flow field, the false diffusion problem in the 
enthalpy advection computation (equations (2.55)) becomes problematic in the 
regenerator and leads to significant errors in the simulated heater and cooler 
heat transfers. This indicates the appropriateness of activating a 
regenerator enthalpy transport model of the kind defined by figure 3.6. In 
this model, linear upwind spatial extrapolations of the volume -aver aged 
temperature field lead to a better approximation of the actual advected 
temperatures within the regenerator than those determined from equations 
(2.55). The particulars of the model are given by equations (3.4) as follows 
(the averaging notation has been dropped for clarity): 
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Table 3.3 Simulation non-dimensional parameters and empirical correlation factors 
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Over the entire working fluid path, not just the heater. 




Figure 3.6 Regenerator enthalpy transport model. 


if (£n>i - 0 then: 

< r n)i - ^(1 + 0.5 Kril/ln)^) - 0.5 K r (l/l n ) (3.4.1) 
if ( g n )± < 0 then: 

(J n ) i - T^l + 0.5 K r li/(l n ) Ul ) - 0.5 * r (Ii/<l n > ln >r 1+1 (3.4.2) 


The upwind extrapolation format of these equations is necessary to maintain 
the transportive property of equations (2.55). The porous advection 
coefficient K r is introduced as a means of accounting for the deviation of 
actual regenerator behavior (for example, non-linear matrix temperature 
profiles and oscillating flow effects) from the ideal behavior suggested. 
Equations (3.4) apply at all the regenerator discrete momentum volume 
centroids (NRF to NRL+1) with the following exceptions: 


if (gn^NRF - 0 then: (T n ) m F - r NHL 

if < 0 then: (^h^nrl+i “ ^nkf 


( 3 . 5 . 1 ) 

( 3 . 5 . 2 ) 


It may be noted that at discrete momentum volume centroid NRF+1 , the adjacent 
heater temperature would be used for T±_ z in equation (3.4.1). Similarly, at 
discrete momentum volume centroid NRL, T ^ is substituted for T^ +1 in equation 
(3.4.2). This is felt to be physically more consistent than the alternative 
of assuming that T NHL occurs at A and T ^ occurs at B in figure 3.6. In 
practice though, the difference between the two approaches is minimal. 
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Owing to the magnitude of the maximum Reynolds numbers occurring over 
the flow field (table 3.3) in the SPDE, it was found necessary to set K r 
almost at unity to obtain agreement with the experimental heat transfer data. 
However, it is important to note in the context of comparing the one- and two- 
dimensional heater module simulations that no modifications are made to the 
nominal Kays and London empirical correlations in the heater. 

The maximum Mach and Reynolds numbers are similar for all three 
simulations. The smallness of the Mach numbers attests to the absence of any 
choking while the magnitude of the Reynolds numbers is indicative of the level 
of flow turbulence. Of particular importance are the low characteristic 
numbers (N ch ) of 24 to 26 (compared, for example, with a N c ^ range of 60 to 
250, depending on operating parameters, for the GM-GPU3 engine (Go87)). This 
indicates that with only about 25 complete information traverses per cycle, 
ignoring information propagation effects in a simulation may not be 
automatically justified. This is borne out by the comparison of the 
experimental and simulation data given in table 3.4 for. the two experimental 
test points. 

Generally, the one -dimensional UPD simulation and experimental results 
are in agreement. The maximum energy balance discrepancy is less than 5.5% 
(external heat supply for test 42) while the simulated and measured expansion 
and compression space mean cyclic temperature and pressure parameters are in 
reasonable agreement. In contrast, the one -dimensional equilibrium simulation 
shows an overall discrepancy of about 30% in the indicated piston work. The 
major source of this error is the mismatch between the measured and simulated 
compression space pressure profile phase angles (which are given relative to 
the piston displacement) . Since the indicated wor.: is proportional to the 

of the phase angle , a 1.2° discrepancy makes a 19% contribution towards 
the indicated piston work discrepancy. It should be noted that no energy 
balance errors are reported for the simulations owing to a data output 
processing error discovered in the code. This error was discovered and 
corrected during the METR simulation runs reported in chapter 5. The 
corrected code typically yields energy balance errors of the order of .01%. 

Hence it is possible that the indicated work discrepancy between the 
equilibrium and UPD algorithms is a consequence of information propagation 
effects. This is suggested by the relative agreement between all the 
remaining experimental and one -dimensional equilibrium simulation parameters 
with the exception of the mean compression space temperature. This is higher 
in the simulation because of the larger predicted pressure profile phase 
angle. However, the overall energy balance simulation results produced by the 
equilibrium algorithm (for test 46 at least) apparently agree with those 
produced by other simulation codes (Te88) such as the Gedeon GLIMPS (Ge86) and 
NASA Lewis SNAP (Te83) codes. In contrast, the MTI harmonic analysis code 
appears to conform to the UPD algorithm predictions (Te88) . However, this 
superficic 1 inter-code comparison is probably only of anecdotal significance 
because, in this context, a detailed irreversibility comparison at least is 
necessary in order to understand how well the codes compare with each other. 
Unfortunately, such data is not yet available. 
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Table 3.4 Experimental/simulation data comparison for the one- and two-dimensional heater modules. 
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Taken to be the mean of the left- and right-hand engine test results. 

Calculated energy balance errors are too large owing to a data output processing error, 
errors produced by corrected code are about .01% at most (see table 5.3). 

Mean of peripheral and central lower compression space values. 



In the face of this dilemma, several other suggestions have been made 
to explain the inaccuracy of the equilibrium algorithm, for example, gas 
leakage between the expansion and compression spaces which is not explicitly 
modelled. In particular, if the inclusion of gas leakage were to eliminate 
the equilibrium algorithm error, then it may be concluded that information 
propagation is not important and that the apparent accuracy of the UPD 
algorithm is a spurious numerical effect. 

In the light of the uncertainty about the information propagation issue 
and the apparent agreement between the equilibrium algorithm, GLIMPS and NASA 
codes for the SPDE, NASA and the principal investigator felt that the two- 
dimensional module simulation should proceed using the equilibrium algorithm 
only pending further work on the information propagation issue, which is 
discussed in chapter 4. 

The predictions of the equilibrium simulations for the one- and two- 
dimensional heaters are similar, at least from an overall cyclic performance 
perspective For test 46, the simulated external heat supplied using two- 
dimensional heater module is 2.6% larger than that predicted using the one- 
dimensional heater module (with the latter value being within .03% of the 
measured external heat supply). In contrast, for test 42, the two-dimensional 
simulated external heat supplied is in closer agreement with the measured 
value than the corresponding one -dimensional prediction. The expansion space 
mean cyclic temperature (which is strongly influenced by the heat supplied in 
the heater) does not exhibit contradictory behavior since, for both tests, the 
two-dimensional simulated value is less than the one -dimensional simulated 
value which in turn is less than the experimental value. 

In this light, a more detailed view of the influence of two-dimensional 
effects may be discussed in terms of figures 3.7 to 3.10. Figures 3.7 and 3.8 
show the cumulative heater wall/fluid heat transfer as a function of angle for 
tests 46 and 42 respectively while figures 3.9 and 3.10 show the volume - 
averaged midpoint heater fluid temperature profiles. In all cases, the one- 
dimensional heater UPD and equilibrium algorithm profiles are in close 
agreement. The two-dimensional heater heat transfer profiles reflect lower 
cumulative heat transfers throughout the cycle, converging rapidly towards and 
then exceeding the one -dimensional profiles beyond 320? The two-dimensional 
volume -averaged heater midpoint temperatures are less than their one- 
dimensional counterparts within a range of about 5 to 10 K over the cycle. 
Hence it is evident that the mechanism of heat transfer simulated in the one- 
dimensional heater using empirical correlations is different from that 
simulated in the two-dimensional heater without such empiricism. In the light 
of the high cyclic Reynolds numbers and resultant turbulence, the numerical 
discrepancy is an expected result since the two-dimensional simulation assumes 
laminar flow in the heater. Thus inclusion of a turbulence model in the 
heater should at least reduce the instantaneous numerical discrepancy, 
although the extent of the narrowing depends not only on the efficacy of 
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the turbulence model but also on the errors inherent in using a steady- state 
correlation. These issues can only be satisfactorily addressed by an 
experimentally based investigation such as that of the METR (see chapter 5) . 
However, it must be noted that from an overall or systems cyclic energy 
balance perspective, the one- and two-dimensional heater module simulations do 
not show any appreciable differences. This emphasizes the limitations of 
inferring the validity of a simulation code purely from time-averaged 
parameter and cyclic energy balance data. Apparent simulation validation 
based on such criteria may be a result of error cancellation effects rather 
than an indication of actual ^e accuracy. 

The aggregate system velocj.cy fields produced for test 42 for the one- 
and two-dimensional heater modules are depicted in figures 3.11 and 3.12. In 
each case, the cyclic angle is plotted on the X axis while the location along 
the one -dimensional working fluid path is denoted by the Y axis. In terms of 
the sign convention of figure 3.1, the expansion space is located on the left 
(Y ss 0) and the lower compression space on the right, the two spaces being 
separated by the heater, regenerator, cooler, upper compression space and 
conical connecting passage. The dashed lines indicate negative velocities. 
Figures 3.11 and 3.12 exhibit no qualitative differences and are similar in 
shape. Quantitatively, there are small differences between the minimum and 
maximum velocities amounting to 3% at most. The test 42 one- and two- 
dimensional heater module simulation aggregated temperature fields shown in 
figures 3.13 and 3.14 exhibit comparative behavior similar to that of the 
velocity fields. The notable, yet consistent, exception is that the maximum 
temperature in the one -dimensional heater module is 9K greater than that in 
the two-dimensional heater module as shown in figure 3.10. 

These aggregate system parameter profiles again show that no major 
qualitative differences are introduced into the simulation by using a two- 
dimensional heater module even on a transient basis when assuming laminar flow 
in the heater. Quantitative effects are also relatively small and limited to 
the heater itself. Therefore, it seems reasonable to suppose that system- 
imposed boundary conditions on the two-dimensional heater module force the 
two-dimensional flow in the module to conform on aggregate to that simulated 
using a one -dimensional heater module. This in turn suggests that two- 
dimensional component simulation may not be an effective way of improving the 
system accuracy of Stirling machine simulation. Such improved accuracy is 
probably only realisable by simulating most (if not all) of the working fluid 
path in two -dimens ions . Details of the two-dimensional flow field in the 
aggregated heater tube are not presented here since they do not bear directly 
upon the systems nature of this discussion. It is sufficient to note that 
these profiles are in agreement with laminar oscillating flow results 
published in the literature. Details of the two-dimensional flow field are 
presented in chapter 3 where they are evaluated against experimental data 
produced by the METR. 
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Figure 3.11 SPDE velocity field (equilibrium algorithm. 1-dim. heater) 
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Figure 3.12 SPDE velocity field (equilibrium algorithm, 2-dim. heater) 
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Figure 3.13 SPDE temperature field (equilibrium algorithm. 1-dim. heater) 
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Figure 3.14 SPDE temperature field (equilibrium algorithm, 2-dim. heater) 









3.6 


CLOSURE 


The baseline comparison of the one- and two-dimensional heater module 
simulated and experimental data shows that two-dimensional flow effects do 
impact the transient heat transfer predictions in the heater itself. However, 
these effects are not proportionately manifest on the system level where t ey 
only minimally affect the cyclic energy balance. Turbulence effects are 
likely to be the major cause of the lack of quantitative correspondence 
between the simulated flows in the one- and two-dimensional heater modules 
since the former nominally includes the impact of turbulence via heat transfer 
coefficient and friction factor correlations. While this justifies isolated 
component studies of oscillating turbulent flow, it does not alter the 
contention that an overall improvement in Stirling machine simulation accuracy 
probably is achievable only if a major portion of the fluid flow path (at 
least the heater, regenerator, and cooler assembly) is simulated m two- 

dimensions . 


The information propagation issue has been established as a possible 
performance issue in SPDE class (high frequency, high pressurisation) Stirling 
engines . This subject is by no means new to Stirling machine simulation an 
has been dealt with by others (such as Organ (Or82) and Taylor (Ta84)) 
although chiefly in the context of method of characteristics simulations e 

information propagation issue has engendered much contention among Stirling 
machine analysts with some suggesting that it is not physically relevant, 
the results produced here serve no other purpose than to stimulate debate and 
definitive research to resolve the information propagation issue, they wi 


have served NASA well . 
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CHAPTER 4 


COMPARISON AGAINST AN ANALYTIC SOLUTION 


4.1 INTRODUCTION 


One of the conclusions from the simulation of the SPDE has been that, 
apparently, an accounting of information propagation effects permits 
convergence between the simulated and measured cyclic energy performance 
parameters to be achieved. In particular, the simulated piston indicated work 
can be made to agree with its measured counterpart only when the Unitary 
Pressure Domain (UPD) algorithm is used. Invocation of the equilibrium 
algorithm consistently leads to an over-prediction of the piston indicated 
work. This over-prediction is principally related to the phase angle between 
the piston displacement and compression space pressure profiles. The 
equilibrium algorithm produces an over-estimate of the phase angle which, 
although small in absolute terms (1 to 2 degrees), is significant because the 
indicated work is proportional to the sine of the phase angle. This is the 
major contributor towards the observed simulation/experimental discrepancies 
of about 30% . 

The following postulates may be used to explain the apparent fidelity 
of the UPD algorithm: 

The information propagation effect is physically significant and 
is being portrayed accurately by the UPD algorithm. 

The information propagation effect is physically significant and 
its portrayal by the UPD algorithm is a spurious numerical effect. 

Information propagation is not significant and the apparent 
accuracy of the UPD algorithm is fortuitous. 

A feasible means of testing these postulates is to apply the UPD and 
equilibrium algorithms to a problem that has a well-defined and validated 
closed- form, analytic solution. Ideally, this problem should closely 
correspond with the boundary conditions prevailing in Stirling machines in 
general and in the SPDE in particular. 

The process adopted has been to apply the existing UPD and equilibrium 
algorithms as used in the one -dimensional SPDE simulation to the selected 
problem. This enables an assessment of whether information propagation 
effects offer a physically reasonable explanation for the observed behavior of 
the algorithms when applied to the SPDE. Thereafter, the numerical accuracy 
of the algorithms is assessed in order to investigate whether spurious 
numerical effects are occurring. 
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4.2 


PROBLEM SELECTION 


A general class of analytically soluble problems involving information 
propagation effects revolves around the description of acoustic transmission 
phenomena in plain tubes. A classic treatise on these problems by Lord 
Rayleigh (St26) encompasses analytic solutions including non-linear effects 
such as viscous dissipation and finite boundary heat transfer. 

Using Rayleigh's development of Kirchoff's equations of sound, Iberall 
(Ib50) developed a first order analytic solution for the transmission line 
problem. This problem is defined geometrically by a length of tube connecting 
an infinite cavity with a rigid receiver volume. The analysis seeks to define 
the phase lag and amplitude ratio of the pressure profile in the receiver 
volume with respect to a sinusoidal pressure variation within the cavi y. 
Physically, the problem is representative of a transmission line connecting a 
signal source to a pressure transducer. The desired outcome of the analysis 
is a means of providing design guidelines on the length and diameter of the 
transmission line so as to ensure adequate measurement accuracy of the 
pressure transducer . 

In another investigation, Chester (Ch64) analyzed the behavior of 
resonant oscillations in closed tubes. The prescribed boundary conditions are 
that one end of the tube is closed while the other is excited by a piston 
oscillating at near- resonant frequencies. The objective of the analysis s o 
investigate the impact of compressive viscosity and boundary shear viscosity 
on the gas oscillations in a well-defined frequency band around resonance in 
which shock waves occur . 

In another analysis, Jimenez (Ji73) extends Chester's analysis to a 
case in which the closed end of the tube is replaced with an arbitrary closure 
condition ranging from fully open to fully closed By expanding the momentum 
equations in terms of a Mach number series, both the amplitude and form of the 
oscillations are predicted in order to show that in both the fully open and 
fully closed cases, shock waves are needed to describe the observed resonant 

behavior. 

Although the problem described by Chester and Jimenez is physically 
more in conformity with the SPDE geometry, their analyses do not correspond 
with the general flow situation prevailing in Stirling machines, particularly 
owing to their appropriate neglect of heat transfer effects . Furthermore 
these analyses are focussed on resonant effects rather than on conditions far 
from resonance where the classical linearised theory is assumed to prevail. 
Hence despite its lack of boundary condition conformity with Stirling 
machines Iberall 's analysis is preferred as a means of investigating the 
validity * of the aforementioned postulates. This is justified by the inclusion 
of a more complete fluid dynamic treatment in the analysis and its yield of 
information more amenable to intuitive interpretation because of the practical 
engineering relevance of the transmission line problem. Furthermore several 
experimental validations of Iberall' s analysis have been performed (Wa65 
G 068 ) which lend credence to using the analysis as a benchmark against which 
the simulation may be compared. 
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4.3 


IBERALL'S ANALYSTS 


The geometry of the transmission line is shown in figure 4.1. 



Figure 4.1 Transmission line geometry. 


A tube of length 1 and diameter d is connected between a cavity and an 
instrument volume . The following assumptions that are significant to the 
comparison exercise are made: 

1. The walls of the instrument volume are flexible so that the actual 
volume may be replaced by a larger equivalent volume that will 
store the same mass of fluid per unit of pressure change. 

2. The flow is laminar throughout the tube. 

3. The gas expands and contracts isothermally in the instrument 
cavity. 

4. The excitation oscillatory pressure is sinusoidal. 

Within these constraint s, the following factors are taken into account in the 
analysis: 

compressible flow in the tube 
finite excitation pressure amplitudes 
fluid acceleration 
finite length of tubing 
boundary heat conduction 

The analytic solution is represented by equation (116) (not repeated 
here) in Ib50. Owing to its complexity, the equation may be conveniently 
solved via a computer program. A key paft of the solution involves the 
computation of Bessel functions. In the computerized implementation of the 
analysis offered by Watts (Wa65), a normal series solution is used for kinetic 
Reynolds (or Valensi) numbers less than 200, while an asymptotic series 
solution is used otherwise. Goldschmied (Go68), however, apparently uses a 
normal series solution throughout. In the implementation used here, Watts' 
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approach is preferred due to its demonstrated superior numerical accuracy. 
However while Watts performs the calculations in real arithmetic, complex 
arithmetic is used throughout here (following Goldschmied) in order to avoid 
what appears to be a coding inconsistency in Watts' implementation. 

In both Watts' and Goldschmied' s implementations, no details of the 
method used to determine the excitation pressure amplitude are given. As the 
amplitude is a critical parameter from a simulation point of view, the 
methodology adopted here has been to choose approximately the largest 
amplitude that produces simulated laminar flow throughout the tube (Reynolds 
numbers less than 2000) . Under certain conditions (particularly at higher 
frequencies) this still yields excitation pressure amplitudes that are very 
small in comparison with the mean pressure. This produces some unavoidable 
truncation errors in the simulation, which although minimized through the use 
of double precision arithmetic, must be borne in mind when interpreting the 
results produced. 


4.4 THE SIMULATION MODEL 

The standard SPDE simulation programs have been applied without 
alteration to the geometry defining a transmission line. The codes have been 
modified to accept a pressure profile and a rigid cavity as boundary 
conditions which replace the piston/cylinder boundaries of the SPDE codes 
The fluid dynamics of the instrument cavity are modelled in full so that the 
isothermal assumption made in Iberall's analysis is not replicated by the 
simulation. The standard Kays and London (KL64) friction factor and heat 
transfer coefficient correlations used in the SPDE simulations are applied 
without any empirical corrections. Care has been taken to ensure that both 
Iberall's analysis and the simulation use identically the same thermodynamic 
initial conditions. This has required some reworking of the constants in 
Iberall's analysis to reflect the use of temperature rather than density as an 
initial condition specifier. 


In view of the varying temporal resolution produced by the UPD and 
equilibrium algorithms, extracting amplitude and phase information from the 
simulated pressure profiles can be subject to large errors, particularly at 
low characteristic numbers (less than about 20). Several approaches have been 
tested including sine transforms, sine quadratures and globally implicit cubic 
spline fits The last approach has produced the most reliable results and has 
therefore been adopted. At each data point, the known excitation amplitude 
and phase angle have been used to check the integrity of the cubic spline 
fits. 


4.5 APPLICATIONS 

Two case studies have been performed to compare the analytical and 
simulation results. The first geometry corresponds to that used by Watts 
(Wa65) to experimentally validate Iberall's analysis. This geometry is 
preferred to that of Goldschmied (Go68) in view of the good agreement of 
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Watts' experimental data with that calculated, especially under conditions of 
resonance . 

The second geometry represents a crude approximation of the SPDE in 
which the flow passage between the expansion and compression spaces is 
modelled as a length of tube with a diameter equal to the average of the 
heater and cooler tube diameters. The end effects of the variable volume 
spaces are neglected so that the tube termination geometry is identical to 
that used in Watts' experiment. 


4-5.1 Watts' Experimental Rig 

The invariant parameters used to describe the geometry and initial 
conditions of the transmission line are listed in table 4.1. These parameters 
correspond with those given in Wa65 to describe the 4- ft transmission line 
tested. The additional input parameters required, namely, the excitation 
pressure amplitude and frequency, are treated as variables. 


Table 4.1 Transmission line parameters. 


Transmission line length (mm) 1220 

Transmission line internal diameter(mm) 7.818 

Instrument cavity diameter (mm) 14 

Instrument cavity volume (cm A 3) .414 

Working fluid ^ir 

Mean excitation pressure (bar) .785 

Mean system temperature (deg C) 25.745 


Using a constant excitation pressure amplitude of 0.016 bars, Iberall's 
analysis yields the results plotted in figures 4.2, 4.3, and 4.4. Figure 4.2 
shows the exponential decrease in characteristic number as a function of 
excitation frequency. Figure 4.3 depicts the variation in amplitude ratio 
(that is, the ratio of the instrument cavity pressure profile amplitude to 
that of the excitation pressure profile) with excitation frequency. Figure 
4.4 shows the corresponding variation in the phase lag of the instrument 
cavity pressure profile to that of the excitation pressure. The amplitude 
ratio profile replicates that reported by Watts and hence shares the 
experimental validity of Watts' data (Watts does not report any phase lag 
comparison data) . The data span an excitation frequency range centered upon 
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FIGURE 4.4 






the first resonant peak. The topologies of figure 4.3 and 4.4 are typical of 
acoustic phenomena occurring in tubes. Of particular note is that the first 
resonance occurs at an excitation frequency of approximately 66.6 Hz, which 
corresponds to a characteristic number of about 4.3 . This may be compared 
against a characteristic number of 2 (Ch64) predicted by acoustic theory. The 
diff ererice in characteristic numbers shows the impact of the irreversibilities 
and non-linearities included in Iberall's analysis. 

Using figures 4.2 and 4.3 as a guide, the simulation has been exercised 
at analytic characteristic numbers of 5.6, 10.0, and 24.1 . These values 
correspond, respectively, to conditions close to resonance, midway between 
resonance and quiescence, and far from resonance. The comparison between the 
simulation and analysis at a characteristic number of 5.6 is summarized by 
figures 4.5 and 4.6. In these figures, the amplitude ratio and phase lag are 
plotted as a function of the number of integration increments per cycle. In 
this context, the lower bound of the variation range corresponds to the number 
of integration increments produced by the UPD algorithm (that is, the first or 
leftmost point plotted is generated by the UPD algorithm) while all the other 
points are generated using the equilibrium algorithm. The demarcated profiles 
are generated by the simulation while the unmarked horizontal lines depict the 
analytical values. 


Figure 4.5 clearly shows that the simulated amplitude ratio asymptotes 
toward about 2.3 with increasing temporal discretisation (or number of 
increments per cycle). At 200 increments/cycle the discrepancy between the 
simulated and analytical amplitude ratio amounts to 14%. However, the UPD 
algorithm (6 increments/cycle) produces an amplitude ratio discrepancy of 
almost 54%. An examination of figure 4.6 shows a similar trend for the phase 
angle behavior in which the simulated values asymptote toward about 5 8 deg 
yieiding a discrepancy of 48.4% with the 11.251 deg phase lag produced by 
Iberall s analysis. The UPD algorithm produces phase lags of 38 and 34 deg at 
t e instrument cavity pressure profile minimum and maximum, respectively 
indicating a mean discrepancy of 220%. 


These comparison data generated under conditions close to resonance 
indicate that two mechanisms are operating simultaneously, one numerical and 
one physical. Numerically, at low characteristic numbers (below 20) the 
numerical accuracy of the UPD algorithm is insufficient to represent’ 
adequately the transient phenomena taking place. In other words, there is a 
certain minimum temporal discretisation below which use of the UPD algorithm 
is not warranted as a result of its numerical accuracy limitations. 
Physically, information propagation effects similar to those described by 
Jimenez (Ji73) are occurring. These effects result in a progressive 
steepening of the pressure wave until such time as at the onset of resonance 
a discontinuity, or shock, forms. Because both the equilibrium and UPD 
information propagation algorithms do not describe shock formation phenomena 

/co^ g0rL are invalid in the presence of pressure wave discontinuities 
(Go87) . Hence, should resonance -induced pressure wave discontinuities exist 
m the SPDE, for example, then neither the UPD nor the equilibrium algorithm 
could, by definition, correctly model the information propagation However 

in view of the SPDE having N ch > 20, the existence of such discontinuities is 
considered unlikely. 
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Figures 4.7 and 4.8 depict the amplitude ratio and phase lag behavior 
at a characteristic number of 10 (excitation amplitude and frequency of .022 
bars and 23.8 Hz respectively). The profiles show the same topologies as 
those of figures 4.5 and 4.6 except that the phase lags at the instrument 
cavity pressure profile maximum and minimum appear to be more divergent 
because of the plot scaling. The asymptotic phase lag discrepancy is reduced 
to about 33% from the 48% noted under conditions closer to resonance. 

Similarly, the asymptotic amplitude ratio discrepancy is reduced to 4%. Hence 
the observations made above are applicable under these conditions as well with 
the additional note that the influence of information propagation effects 
becomes less pronounced the further the operating point from resonance. 

Lastly, at a characteristic number of 24.1 (excitation amplitude and 
frequency of ’.05 bars and 11.7 Hz, respectively), figures 4.9 and 4.10 show 
that conditions are sufficiently far from resonance that only numerical 
effects remain. In figure 4.9, the amplitude ratio discrepancy varies between 
1.1% and .9% for the UPD and equilibrium algorithms, respectively. Figure 
4.10 effectively shows that the simulated phase lags bound the analytic phase 
lag in a somewhat random fashion which, owing to their smallness, is 
symptomatic of numerical effects. As the SPDE also yields a characteristic 
number of 24, these data show that the temporal accuracy of the UPD algorithm 
is just adequate to portray the information propagation behavior of the SPDE 
with some qualitative confidence although, quantitatively, the 
representational accuracy of the results is probably somewhat less than that 
suggested by the small SPDE simulation energy balance errors achieved. 


4.5.2 Pseudo -SPDE Geometry 

The invariant parameters used to define the pseudo-SPDE geometry are 
listed in table 4.2. 


Table 4.2 Pseudo-SPDE parameters. 


Transmission line length (mm) 
Transmission line internal diameter(mm) 
Instrument cavity diameter (mm) 
Instrument cavity volume (cm A 3) 

Working fluid 

Mean excitation pressure (bar) 

Mean system temperature (deg C) 


321.25 

1.397 

14 

.414 

Helium 

150.29 

100 
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In selecting this geometrical configuration, no attempt has been made to take 
account of the effect of either the expansion or the compression space by 
increasing the volume of the instrument cavity. Although Iberall and 
Goldschmied report data for volume ratios of as much as 50 to 1, in the case 
of the SPDE expansion space for example, using the tube data of table 2, the 
cavity- to -tube volume ratio amounts to about 522 to 1. In terms of the 
validation data reported by Watts and Goldschmied as well as the assumptions 
made in Iberall's analysis (see section 4.3), there does not appear to be any 
solid basis for extrapolating the analysis by an order of magnitude in volume 
ratio terms. Furthermore, the mechanism of pressure wave generation in the 
SPDE is very different from that used by Iberall. Hence the main value of 
Iberall's analysis in terms of understanding the gas dynamics of the SPDE is 
to test whether closed tube resonance effects are a likely cause of the 
observed phenomena. In this context, it thus seems desirable to maintain the 
congruency of the tube boundary conditions with those that produced the high 
level of agreement between Watts' experimental data and Iberall's analysis. 

Using an excitation amplitude of .9 bars, the results of applying 
Iberall's analysis to the pseudo -SPDE geometry are shown in figures 4.11 to 
4.13. Figures 4.11 and 4.12 show that the first resonant peak occurs at an 
excitation frequency of 504.8 Hz, 5.1 times larger than the nominal 100 Hz 
operating frequency. At 100 Hz, the amplitude ratio amounts to 1.05 while the 
phase lag is .327 deg (figure 4.13) and the characteristic number is 35.2 . 
Hence even though the pseudo- SPDE geometry characteristic number is greater 
than that of the actual SPDE (approximately 24) , these results show that 
closed tube resonant effects are not likely to be the cause of the observed 
SPDE phase angle behavior (the difference in characteristic numbers serves to 
increase the conservative nature of this inference) . 

A comparison of the simulated and analytical data for the pseudo-SPDE 
geometry is given in figures 4.14 and 4.15. The comparison is performed at a 
frequency of 99.569 Hz (taken from an MTI test) and an excitation amplitude of 
.9 bars. Figure 4.14 yields an amplitude ratio discrepancy of .55% for the 
UPD algorithm, increasing to .65% for the equilibrium algorithm at 90 
increments per cycle. It is interesting to note that in this case the 
simulation over-predicts the amplitude ratio while in the previous case (in 
which the transmission line is 3.8 times longer) the simulation under-predicts 
the amplitude ratio. This is a probable result of the neglect of the 
advection terms in the momentum equation which Iberall uses to produce 
equation (116) (Ib50). V 

[In the second paragraph of section 4 on page 100 of Ib50, Iberall makes 
reference to the fact that his analysis is based on Kirchoff's equations of 
sound. Tracing this to the relevant equation quoted by Rayleigh (St26 paee 
315, equation 13), in the class of problems being considered, Iberall 
apparently agrees with Rayleigh's assertion (St26, page 3) that Whenever 
the motion is very small, the (advective) terms u du/dx, etc., diminish in 
relative importance, and ultimately, D/Dt - d/dt (or the substantive and total 
envatives are equal)." This assumption does not apply to SPDE conditions 
nor does it lend credence to cases in which large instrument cavities produce 
rluxes big enough to produce significant momentum transport. These 
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limitations should be borne in mind when comparing the simulated results 
(which include all the advective terms) with their analytic counterparts.] 

Due to the smallness of the phase lags, figure 4.15 is inconclusive in 
revealing any trends that cannot be ascribed to systematic numerical effects , 
in particular those produced by the cubic spline fitting procedure. 


4.6 INTERPRETATION 


Under resonant conditions in which information propagation effects are 
responsible for the observed transmission line behavior, the inadequacy of the 
equilibrium information propagation hypothesis is demonstrated. This is 
manifested by the equilibrium algorithm producing results that asymptote 
towards values which are not in agreement with those produced by Iberall's 
analysis. Hence, should such effects exist in the SPDE, they would not be 
properly accounted for by an equilibrium information propagation hypothesis. 

As both the UPD and equilibrium algorithms are based on a temporal 
integration scheme which is only first order accurate, the comparison data 
show that the algorithms quickly lose their physical validity at temporal 
resolutions below 20 to 30 increments per cycle. Hence, use of the UPD 
algorithm to simulate the SPDE at an effective temporal resolution of 23 to 24 
increments/cycle is somewhat risky as the integration algorithm is operating 
in a region of marginal temporal accuracy. Therefore, the comparison exercise 
is inconclusive in determining whether the UPD algorithm accurately can 

portray the information propagation behavior of the SPDE in a quantitative 
sense . 


In other words , the data show that the numerical errors produced by a 
first order accurate, implicit integration algorithm at time steps large 
enough to satisfy the UPD hypothesis at low characteristic numbers do not 
permit an assessment of the ph ysical accuracy of the UPD algorithm to be made 
This conclusion in effect reflects the limitations of the pressure- linking 
algorithm itself, which is defined in essence by equation (2.47). Since the 
off-diagonal pressure term coefficients are given as a product of the time 
step At and a geometrical factor A z /V, these coefficients are positive 
definite and are independent of the working fluid state. Therefore, an 
implicit solution of the pressure equation (2.47) implies infinite information 
propagation irrespective of the time step. Thus, satisfying the UPD 
hypothesis is necessary physically in order to ensure that the actual and 
numerical information propagation rates match. Decreasing the time step 
towards the equilibrium limit thus effectively increases the information 
propagation rate towards infinity. Therefore the data show that the numerical 
effect of the time step size predominates, submerging the simultaneously 
occurring physical effect. 

Thus before any further assessment of the physical relevance of the 
information propagation hypothesis is attempted, the impact of numerical 
accuracy on the UPD algorithm in particular needs investigation. 
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4.7 


NUMERICAL ACCURACY CONSIDERATIONS 


The physical adequacy of the first order temporal accuracy of the 
simulation code at low characteristic numbers may be evaluated by using a more 
accurate temporal integration scheme. The issue of significance in this 
endeavour is to find a temporally more accurate scheme which does not destroy 
the hypothesized information propagation modelling characteristics of the U 
algorithm in particular. After studying the Gedeon GLIMPS methodology (Ge86) 
and further review of the literature, the following observations can be made. 


The GLIMPS code is a particular case of a general class of 
numerical algorithms which may be conveniently referred to as 'multi- 
step implicit' methods. In these methods, a given time increment is 
resolved into m sub-steps so that each sub-step is solved implicitly 
via a single matrix inversion, with both m and the magnitude of the 
time increment being arbitrary. Each sub-step m may be further divided 
into n increments so that the order of the temporal gradient 
approximation used in each sub- step is restricted only by the 
requirement that k < n, where k is the number of increments required to 
implement a particular integration scheme. For example, a 5th order 
backward difference approximation would require n > k - 6. 


Hence the GLIMPS code uses a time increment equal to the cycle 
period with one sub-step (that is, m - 1) with n apparently being in 
the range of 6 to 10 and k - 3. In contrast, the 'UPD simulation' has 
n _ _ \ with m set equal to the integer nearest the characteristic 

number. Hence, if the numerical accuracy of multi-step implicit 
methods is a monotonically increasing function of n, then a generalised 
version of the GLIMPS - code (without the cyclic equilibrium closure 
restriction, in particular) theoretically should be physically valid 
for characteristic numbers > 1. However, as a general proposition, 
multi-step implicit methods should be valid for all characteristic 
numbers, although it is hard to conceive of Stirling machines being 
built with characteristic numbers approaching 1. 


b It is possible to construct a fully implicit temporal integration 

procedure based on the Crank-Nicholson scheme, which is second order 
accurate, while retaining the information propagation characteristics 
of the UPD algorithm. 

Another approach to testing the numerical accuracy of the algorithms is 
to implement the full pressure domain splitting (PDS) algorithm discussed in 
section 2.6. However, in view of the large increment in computational effort 
associated with the PDS algorithm and, therefore, NASA's reluctance to pursue 
this approach, it was decided to limit the numerical accuracy investigation of 
the 'UPD simulation' to an implementation of the Crank-Nicholson integration 

scheme . 


4.8 THE CRANK -NICHOLSON INTEGRATION SCHEME 

A generalised two-step temporal integration scheme may be defined by: 
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(4.1) 


<Wd t - K d^/dt (t+At) + (1-JT) d^/dt (t) 

If /C - 0, equation (4.1) represents an explicit scheme of 0(At) ; if K - 
1, equation (4.1) is an implicit scheme of 0(At) ; and, if K - .5 then equation 
(4.1) represents a Crank-Nicholson scheme of 0 (At 2 ) . 

Thus, choosing K - .5 and discretising: 

^(t+At)/4t - .5 d^>/dt( t+ A t j ~ ^(£)/A t + .5 d^/dt(£j (4.2) 

Hence by setting d^/dt^j of the current time step equal to d\i/dt ( £ +A £) 
of the previous time step (which by definition is implicitly determined) , the 
Crank-Nicholson scheme retains much of the stability of the fully implicit 
scheme as confirmed by von Neumann stability analysis (Ro82) , although the 
method is not unconditionally stable. Stability problems arise when At 
becomes large enough for some Fourier modes to overshoot, in which case, the 
0(At ) error terms for the Crank-Nicholson scheme will exceed the 0(At) error 
terms for the implicit scheme. Furthermore, equation (4.2) retains the 
information propagation characteristics of the first order implicit scheme 
used in the UPD simulation. 


4.9 RESULTS 


The parameters describing the transmission line used as the basis for 
the comparison exercise are listed in table 4.3. These parameters correspond 
to those given in Go 68 table V and Wa65 figure 6 . 


Table 4.3 Transmission line parameters. 


Transmission line length (mm) 1220 

Transmission line internal diameter (mm) 4.66 
Instrument cavity length (mm) \ 

Instrument cavity volume (cm A 3) .414 

^ as Air 

Mean excitation pressure (bar) 77 

Excitation pressure amplitude (bar) .00077 

Mean system temperature (deg C) 28.98 


The computerised implementation of Iberall's analysis has been 
carefully rechecked against Watts' results given in Wa65 table 6 to ensure 
exact agreement. 


The simulation parameters used are given in table 4.4 
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Table 4.4 Simulation parameters. 


No. of transmission line contr( volumes 
Global friction factor corre 1 ltiplier 

Global heat transfer coefficit u lation multiplier 

Convergence energy balance error (.%) 


20 

1 

1 

<.05 


The results of exercising the UPD algorithm using both the implicit and 

since N ch is inversely proportional to frequency) . The tabulated data are 
plotted in figures 4.16 to 4.19. 

In all cases, the amplitude ratio errors produced by the Crank- 
Nicholson scheme (CN-s) are less than those produced by the implicit scheme 
! phase angle errors are also less than their I -s counterparts 

for N h > 8.15. Below this, the time steps are large enough so that the CN- 
truncation errors exceed those of the I-s as noted in section 4.8 above. 

It is also evident that the phase angle differences between the C N-« 
and I-s are larger than the amplitude ratio differences, which is consistent 
with the notion that information propagation effects are more apparen , 
numerically, in the phase angles. 

If one accepts that a 10% discrepancy represents a reasonable level of 
validation (see Go87, section 7.7.4), then based on the errors reputed in 
table 45, the CN-s is valid for N ch > 13.27 while the I-s is valid for N ch - 
20 38 . The accuracy differences are also far more pronounced for the phase 
angles than for the amplitude ratios as shown by figures 4.17 and . . 

Based on this discussion, the UPD and equilibrium algorithms are 
compared at N ch of 8.15, 13.27, and 20.38, respectively corresponding to cases 
in which neither the CN-s nor the I-s schemes are valid, only the CN s 
ialld and both schemes are valid. The results are given In tables 4.6 and 
plotted in figures 4.20 through 4.31. 

In all cases, amplitude ratio errors decrease with decreasing 
integration time step (increasing number of increments per cycle). However, 
thfampUtude ratios asymptote to a constant value beyond 50 increments per 
Self w“h the exception of figures 4.29 and 4.31, the phase angle errors 
also reveal this asymptotic behavior except that the lining etrorTs about 
50% for N ch <20. The behavior of the CN-s in figures 4.29 and 4.31 may be 
inconsistent and requires further study before comment. 
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UNITARY PRESSURE DOMAIN ALGORITHM 

Amplitude Ratio Comparison 
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UNITARY PRESSURE DOMAIN ALGORITHM 

Phase Angle Errors 
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Characteristic Number 



Table 4.5 UPD algorithm: Comparison of Crank-Nicholson and implicit integration routines 
Data: Watts figure 6; Goldschmeid table V 
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Table 4.6 Equilibrium algorithm: Comparison 
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No. of Increments per Cycle 


Figure 4.21 
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EQUILIBRIUM ALGORITHM 

Amplitude Ratio Errors at 35 Hz 
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Figure 4.25 
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Amplitude Ratio Errors at 215 Hz 
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Figure 4.27 
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Amplitude Ratio Comparison at 14 Hz 
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Phase Angle Comparison at 14 Hz 
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EQUILIBRIUM ALGORITHM 

Phase Angle Errors at 14 Hz 



Figure 4.31 


In general terms, these phenomena confirm those observed in section 
4.5, namely, that under assumed equilibrium information propagation conditions 
at low characteristic numbers, the errors in modelling information propagation 
using the equilibrium hypothesis produce incorrectly simulated phase angles 
while reducing the simulated amplitude ratio errors. The fact that this 
behavior is exhibited by both the CN-s and I-s indicates the existence of a 
physical effect that is not properly modelled by either the CN-s or the I-s 
techniques owing to their numerical limitations. 


4.10 CONCLUSION 


Essentially it seems that numerical integration errors are a reasonable 
hypothesis for describing the limited range of applicability of the one-step 
implicit application of the UPD algorithm. This conclusion is also consistent 
with the apparent success of the GLIMPS code in duplicating Watts' amplitude 
ratio data for N ch < 8.15 . However, a generalised implementation of a multi- 
step implicit method would have to be developed and applied to the 
transmission line to give this assertion better credence. 

A further complicating issue arises from the difference between the 
UPD/equilibrium and GLIMPS codes in the numerical implementation of the 
conservation balances. The former uses a pressure- linked algorithm (with the 
limitations discussed in section 4.6) while the latter does not. This poses 
the question as to whether a pressure- linked algorithm is valid at low 
characteristic numbers or, for that matter, whether it is valid for gas- 
dynamic flows at all. Hence, there is a possibility that such a lack of 
validity may synergistically or in isolation account for the behavior of the 
UPD/equilibrium simulations when applied to the SPDE as a particular case. In 
general , such behavior has not been observed (Go87) when, apparently, the 
characteristic numbers have been large enough (greater than 60) to make the 
assumption of infinite information propagation reasonable. 

As a practical consequence of this work, it is recommended that the 
SPDE one -dimensional simulation code be converted from an I-s to a CN-s. This 
should provide a better margin of safety between the validity limit of the CN- 
s and the operating characteristic number of the SPDE. 

In terms of the two-dimensional SPDE code, it may be argued that even 
the CN-s UPD simulation is inadequate, particularly in terms of implementing 
workable turbulence models. This arises because of the apparent conflict 
between the smaller time steps required for these models and the larger time 
steps required for modelling information propagation using a first order 
implicit algorithm. Thus, absent the availability of viable alternatives 
(such as PDS or multi-step implicit methods) it seems appropriate at this 
stage to concentrate on the turbulence modelling in the two-dimensional code 
using the equilibrium hypothesis and simply acknowledge the resulting 
existence of information propagation errors at low characteristic numbers. 
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4.11 PROGNOSIS 


In view of the wealth of information gleaned from comparing the 
simulation codes against Iberall's analysis, it is recommended that such a 
comparison be adopted as a universal validity check for Stirling machine 
simulation codes. 

The development of either a PDS or a generalised multi-step implicit 
method (with the latter alternative appearing preferable at this stage) for 
Stirling machine simulation as well as an independent investigation of the 
validity of pressure -linking both appear to be worthwhile future research 
activities , 
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CHAPTER 5 


SIMULATION OF THE MECHANICAL 
ENGINEERING TEST RIG 


5.1 INTRODUCTION 


The philosophy guiding the simulation of the mechanical engineering test 
rig (METR) was to provide a means of using its empirical data in a generalised 
simulation of Stirling cycle machines in both one and two dimensions. The 
following methodology was conceptualised for accomplishing this goal: 

a. Validate a system-based simulation of the METR including a two- 
dimensional description of the test section. The primary goal of 
this exercise was to define a turbulence model that will enable 
the experimental data to be matched in the context of a system 
simulation that includes an accurate description of actual 
oscillating flow boundary conditions. 

b. Based on the validated two-dimensional code predictions, develop a 
set of friction factor and heat transfer correlations which may be 
used in a one -dimensional simulation. 

c. Compare the one -dimensional predictions directly against the 
empirical data and, if necessary, refine the correlations. 

d. Test the one- and two-dimensional simulations against available 
experimental data for actual Stirling hardware, in this case, the 
SPDE . 


The work described in this chapter represents the progress made in 
accomplishing these tasks. The simulation activities listed fundamentally 
depend on the availability of METR experimental data of sufficiently high 
-quality. For this reason, at the termination of the project, sufficient 
experimental data were available only to make a preliminary attempt at 
completing task a (see chapter 1). 

Two sets of experimental data are available. The first set is used to 
define a baseline simulation case against which the effect of turbulence 
models may be measured. This baseline excludes all turbulence modelling from 
the two-dimensional components (that is, the test section is described two- 
dimensionally in purely laminar terms) and uses standard Kays and London 
(KL64) steady- state friction factor and heat transfer correlations for all the 
one -dimensional components. 
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A further objective of the baseline simulations was to provide feedback 
to Simon and Seume on refinements to the METR instrumentation system and, in 
particular, on which additional measurements were necessary in order to 
achieve adequate simulation validation. In this context, it must be noted 
that the baseline data provided by Seume and Simon are of a preliminary nature 
and were intended not only to 'shake down' the METR instrumentation system but 
also to refine the mechanics of the raw experimental data transfer and 
processing. These data are acknowledged to be unsuitable for turbulence model 
development and simulation validation (Se88) . 

The second set of data, which is nominally of sufficient quality for 
simulation development and validation, consists only of half a cycle of data 
at one probe location. Unfortunately, this data set is really too small and 
insufficiently free of systematic errors to allow a truly meaningful attempt 
at validating the turbulence model developed. Hence the model proposed is of 
a tentative nature only and, therefore, fulfillment of the overall goal as 
well as the stated objective (of completely validating the turbulence model 
developed) must await the availability of further experimental data. 


5.2 DESCRIPTION OF THE TEST RIG 

A schematic of the METR as simulated is shown in figure 5.1. The rig 
consists of a combination of seven components arranged axially in a 
symmetrical pattern about the test section midpoint. A piston/cylinder is 
located on the left-hand end of the rig while the right-hand end is open to 
the atmosphere. The universal sign convention adopted is that the positive 
direction is toward the right-hand or open end of the rig. The METR is 
simulated as operating with purely isothermal boundaries since no heat 
transfer measurements are as yet available. Hence all the boundary 
temperatures are kept constant at the ambient value. 

The individual rig component descriptions used in the baseline 
simulations correspond to those listed in table 5.1. This table is the output 
of a geometry definition computer program which enables the rig dimensions to 
be arbitrarily altered in terms of the defined rig components and component 
sequence. The program may be readily modified to include new component 
definitions as well as different sequences of defined components. Table 5.1 
reveals that care has been taken to define those minor details of the rig 
geometry that are judged to have a more than negligible effect on the fluid 
flow. 


The spatial discretisations employed (table 5.1) conform to those 
established as suitable for an integral description of Stirling machine fluid 
dynamics (Go87) . The 11 x 8 mesh used in the test section for the two- 
dimensional description (21 axial control volumes are used in the one- 
dimensional description) is a compromise between adequate spatial resolution 
and the available computing power. It should be noted that 'adequate spatial 
resolution for an integral (or discrete volume) analysis is very different 
from that required for a differential analysis. In this context, increasing 
the spatial resolution toward differential levels would require the use of a 
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Figure 5.1 Schematic of Mechanical Engineering Test Rig 











Table 5.1 METR baseline simulation input parameter set. 


Expansion Space 

Piston stroke (mm) 

Piston end cap thickness (mm) 

Piston end cap taper (deg) 

Clearance length at piston TDC(mm) 
Cylinder diameter (mm) 

Conical Transition Piece 
Axial length (mm) 

Cylinder mating diameter (mm) 
Flexible joint mating diameter (mm) 
__#_nf control volumes 
Flexible Joint 

Cylinder side axial length(mm) 
Cylinder side diameter (mm) 

Cylinder side # of control volumes 
Nozzle side axial length(mm) 

Nozzle side diameter (mm) 

Nozzle side # of control volumes 


355.6 

6.0 

30.0 

3.1750 

355.6 

11.25 

212.725 

355.6 

127.0 

3 


21.0 

152.0 

2 

63.0 

139.7 

2 


Spacer 

Axial length(mm) : 

Cylinder side flange length(mm) : 

Nozzle side flange length (mm) : 

Diameter(mm) • 

Cylinder side flange diameter(mm) : 

Nozzle side flange diameter (mm) 

# of control volumes i 

Heat Exchanger 

Axial length (mm) : 

Cylinder side flange length(mm) : 

Nozzle side flange length(mm) : 

Casing diameter(mm) • 

Hydraulic diameter(mm) : 

Flow area ratio : 

# of control vol umes • 

Nozzle 

Size (1— small , 2-large) : 

# o f half-nozz le control vols : 

Test Section 

External diameter (mm) • 

Axial length(mm) : 

# of axial control volumes (4 minimum) : 

# of radial control volumes (2 minimum) : 

Wall conductivity (W/m-K) : 

Wall density(kg/m A 3) : 

Wall speci fic heat capacitv(J/kg- K) : 

Working Fluid 

__jjorlcing_^luid_Q^l£ A _2^Helium dk _3^l3£dro£en^ i 


63.5 

3.175 

6.35 

139.7 

127.0 

127.0 

3 


50.8 

9.525 

12.7 

127.0 

2.1824 

.556 

4 


1 

3 


50.8 
2609.1 
21 (1-d); 
11 (2-d) 

8 

1.0211 

2466.84 

753.624 

1 
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supercomputer, for which no funds were available. Hence, determining the 
impact of radically changing the test section spatial resolution must of 
necessity be neglected here. 


5.3 BASELINE CASE SIMULATION APPLICATION 


The one- and two-dimensional simulation codes used to model the METR are 
in essence identical to those used for the SPDE. The codes are rearranged to 
conform to the geometry of the METR by excluding the regenerator and the 
compression space and by reconfiguring the connecting ductwork. The two- 
dimensional heater module is adjusted to represent the geometry of the single 
tube test section. The pseudo- two-dimensional heater boundary condition used 
at the expansion space side of the SPDE heater is replaced with a one- 
dimensional boundary condition appropriate to a nozzle (with an axially 
varying flow area) . 

In converting the SPDE codes to the METR geometry, the following 
revisions to the numerical algorithm have been implemented: 

The energy balance reporting error has been corrected yielding 
typical energy balance percentage errors of order 10' 2 or less 
(compared with errors of order 10° previously) . 

Computationally more efficient discretisation methods have been 
developed for calculating the dissipation and reversible 
thermal/kinetic energy conversion terms. 

Provision has been made for describing control volumes that may be 
generated as bodies of revolution based on continuous functions 
(cones, paraboloids, etc.). 

In order to fulfill the objective of generating a baseline laminar flow 
case for the two-dimensional simulation of the test section, the following 
adjustments to the code have been made: 

The Reynolds stress tensor, turbulent heat flux and dissipation 
terms have been removed from the momentum and thermal energy 
equations, respectively (equations 3.3). 

Mass flux fluctuations generated by the iterative integration 
algorithm are numerically filtered out (such fluctuations are a 
manifestation of a chaos model implicitly embedded in the 
algorithm) . A radially polarised cut-off amplitude of 10’ 1 
kg/m 2 .s Is imposed, which has the effect of removing all the 
turbulence with RMS amplitudes larger than the cut-off amplitude; 
the smaller scales remain. 
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5.4 BASELINE CASE RESULTS 

The primary focus of the following discussion is directed towards a 
comparison of the available experimental data with its simulated counterparts. 
This does not constitute a rigorous validation but does serve to expose those 
areas requiring further experimental and analytic attention. 


5.4.1 Experimental Data 

The experimental data set made available by Seume and Simon for 
validating the simulation consists of: 

bulk mean velocities, that is, area-weighted, one -dimensional flow 
velocities 

axial ensemble -averaged mean velocities and root -mean- square (RMS) 
velocity fluctuations at seven locations on two radial planes (at 
the midpoint and toward the cylinder end of the test section) 

These data were gleaned for the parametric conditions shown in table 5.2. 

Table 5.2 Parameters for the baseline case experimental data. 


Date 

Engine Speed 
(rpm) 

Probe Position (mm) 
(relative to cylinder end) 

4/22/88 

11.3 

1304 (midpoint) 

4/24/88 

11.3 

609 

4/26/88 

33.5 

1304 (midpoint) 


The data set supplied is not entirely adequate for simulation validation 
purposes since several key data are missing. Specifically, the missing data 
include : 


ambient pressure and temperature conditions 
transient leakage flow rates past the piston 

The actual rig geometry is considered to be adequately represented by 
the description encompassed within table 5.1 with the usual uncertainties 
present in the fin surface geometry assumed for the heat exchangers. On the 
advice of Seume, a standard Kays and London (KL64) louvered plate -fin surface 
designated ' l/4(b) -11 . 1' has been assumed together with its steady-state 
friction factor and heat transfer coefficient correlations. 

In the absence of any cylinder leakage data, the simulation is performed 
on the basis that no leakage occurs in the cylinder. This may produce 
discrepancies between the experimental and simulated bulk flow velocities. 
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Should cylinder leakage actually occur, then the discrepancies may be used as 
an a posteriori means of assessing what leakage should be incorporated into 
the simulation to match the experimental data. 

The unavailability of measured ambient temperatures and pressures is of 
much greater significance because the rig operates as an open system. 
Furthermore, since the two tests at a constant rig speed of 11.3 rpm were 
performed on different days, local weather conditions are such that 7 to 10 % 
variations in barometric pressure are possible between tests. This can 
produce significant discrepancies in comparing the simulation and experimental 
data. Nevertheless, for the purpose of carrying out the validation exercise, 
it is necessary to assume an ambient temperature and pressure. Values of 
26°C and 0.98274 bar, respectively, were chosen, these being typical averages 
for the days during which the tests were conducted. 


5.4.2 Comparison of System Results at 11.3 rpm 

Figure 5.2 shows a comparison of the experimental bulk velocity profiles 
at 11.3 rpm with those simulated in one dimension. The difference between the 
profiles simulated at probe positions of 609 and 1304 mm are so small relative 
to the scale of the graph that the two profiles appear coincident. Figure 5.3 
compares the experimental bulk velocity at 11.3 rpm at the 1304 mm probe 
position with that produced by the two-dimensional simulation. There is no 
discernable difference between the one- and two-dimensional bulk velocity 
profiles at the 1304 mm probe position, which is also true for the system 
axial velocity profiles shown in figures 5.4 and 5.5. Other than the 
difference in axial discretisation, the one- and two-dimensional simulated 
profiles are identical. 

The phase of the simulated and measured profiles in figures 5.2 and 5.3 
are in agreement, while the difference in amplitude of the profiles shows the 
extent of the piston leakage. This leakage is less than that expected by 
Seume and Simon, who used an ideal velocity profile to make their 
determination. The irregularity of the measured profiles at peak flow 
velocities is believed to be caused by adhesion of the piston seal to the 
cylinder wall and/or other mechanical effects. 

Of particular concern in figure 5.2, however, is the large difference in 
experimental amplitudes relative to the coincidence of those measured. When 
the flow is negative between 0° and 180° (by convention, flow toward the 
piston) , the pressure at 1304 mm is greater than that at 609 mm since the 
pressure gradient is opposed to the shear forces acting on the flow. This is 
confirmed by the both the one- and two-dimensional simulated pressure profiles 
shown in figures 5.6 and 5.7. Hence, since the system is essentially 
isothermal over the test section (see figures 5.8 and 5.9), the density at 
1304 mm is greater than that at 609 mm. Hence, between about 80 and 100 
degrees of crank angle when the flow acceleration is small or zero, in order 
to conserve mass, the velocity at 609 mm should be greater than that at 1304 
mm, which is indeed the case (also shown by the simulated data) . 


103 



SIMULATED AND EXPERIMENTAL BULK VELOCITY PROFILES 

RIG SPEED: 11.3 rpm 
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SIMULATED AND EXPERIMENTAL BULK VELOCITY PROFILES 

RIG SPEED: 11.3 rpm 
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FIGURE 5.4 SYSTEM AXIAL VELOCITY PROFILE FROM 1-D CODE: 11.3rpm 
















FIGURE 5.5 SYSTEM AXIAL VELOCITY PROFILE FROM 














FIGURE 5.6 SYSTEM PRESSURE PROFILE FROM 1-0 CODE: 11.3npm 
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FIGURE 5.8 SYSTEM TEMPERATURE PROFILE FROM 1-D CODE: 11.3rpm 
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FIGURE 5.9 SYSTEM TEMPERATURE PROFILE FROM 2-D CODE: 







Conversely, during positive flow (between 180° and 360°) the velocity at 
609 mm should be less than that at 1304 mm (shown by the simulated data) 
which does not agree with the experimental observations. Furthermore the 
magnitude of the experimental velocity drop is not intuitively credible in 
view of the small pressure drops expected in an atmospheric rig with pisto 
leakage. After consultation with Seume, it was concluded that the most likely 
cause of the discrepancy is a systematic error in the hot-wire anemometer 
calibration. The 609 and 1304 mm tests were performed on different days and 
the anemometer calibration procedures used did not account for barometric 
pressure differences. Hence, because hot-wire anemometry is dependent on 
velocity as well as density, an ambient pressure variation of about 5 to 7 % 
could account reasonably for the apparent axial velocity gradient measured. 

An additional source of systematic error is the manner in which the radial 
mean velocities measured are spatially integrated to yield the bulk flow 
velocity It was thus agreed that future tests should include anemometer 
calibration procedures that include the effect of ambient pressure. Also, a 
better method of computing the bulk velocity was suggested. 

A comparison of the system pressure profiles generated by the one- and 
two-dimensional codes at 11.3 rpm is shown in figures 5.6 and 5.7 The 
pressure drop across the test section predicted by the one -dimensional code is 
larger than that predicted by the two-dimensional code over the cycle. Hence 
the pseudo -steady- state friction factor assumption used in the one -dimensional 
code generates larger dissipation than that produced by the two-dimensional 
code. The differences in profile topology apparent principally at the 
nozzle/test section interfaces are a manifestation of the difference in bulk 
fluid temperatures produced by the one- and two-dimensional codes. The two- 
dimensional code produces a lower bulk temperature over the test section than 
its one -dimensional counterpart (see figures 5.8 and 5.9). This is 
intuitively reasonable since the one -dimensional code, which uses a pseudo- 
steady-state heat transfer correlation, does not describe radial temperature 
gradients. Such gradients are shown clearly by the radial temperature prof e 
(figure 5.10) produced by the two-dimensional code. Over most of the c ycle. 
the temperatures adjacent to the wall are relatively larger than those the 
flow core even though the net temperature difference of .3 K is small. Tims, 
in the two-dimensional code, the actual wall/fluid heat transfer occurs in the 
boundary layer with the radial temperature gradient being constrained chiefly 
by turbulent enthalpy transport. Hence the extent to which the bulk 
temperatures predicted by the one- and two-dimensional codes differ is 
probably a direct effect of the exclusion of such transport from the two- 
dimensional simulation. 


Another difference between figures 5.6 and 5.7 is the relative 
smoothness of the one -dimensional plot. The jaggedness produced by^the two- 
dimensional code is caused by residual simulated turbulence below the cut-off 
of the numerical filter used to exclude such turbulence deliberately. It is 
interesting to note that the occurrence of the low level noise passing the 
filter is approximately in phase with the corresponding high levels of 
measured RMS velocity fluctuation shown in figures 5.11 and 5.12. The 
magnitudes of the pressure drops simulated across the test section are less 
than 500 Pa. An assessment of the accuracy of these predictions requires the 
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FIGURE 5.10 SIMUL. RADIAL TEMPERATURE PROFILE AT 1304mm: U.3rpm 




FIGURE 5.11 EXP. RADIAL PROFILE OF AXIAL RMS VEL. FLUCT. AT 609mm: 11.3rpm 



FIGURE 5.12 EXP. RADIAL PROFILE OF AXIAL RMS VEL . FLUCT. AT 1304mm: 11.3rpm 








installation of appropriate pressure transducers on the test rig. Hopefully, 
this will be accomplished in the future. 

The system temperature profiles of figures 5.8 and 5.9 share a common 
topology but the two-dimensional profile magnifies the deviations from the 
isothermal mean over the test section. The adiabatic temperature variation in 
the cylinder is readily apparent, although the physical significance is 
exaggerated by scale effects (1 to 4 K in 300 K) . 

The measured root-mean-square (RMS) velocity fluctuation profiles shown 
in figures 5.11 and 5.12 reveal the structure of the turbulence occurring in 
the test section. From figure 5.11, it may be observed that when the flow is 
negative, the centerline turbulence is considerably greater than that 
occurring when the flow is positive. This shows the extent to which 
turbulence is advected down the tube in the flow direction. In the positive 
flow direction, the 609 mm probe location is towards the tube entrance where 
the turbulence is nominally small by design. Hence the centerline fluctuation 
is smaller than that occurring during negative flow when the measurement 
location is toward the test section exit. Concomitantly, figure 5.12 shows 
that the centerline velocity fluctuation is largely independent of flow 
direction since the measurements in this case are made at the test section 
midpoint; the advected turbulence is thus the same regardless of the flow 
direction. It is also noteworthy that the largest fluctuations occur close to 
the boundary layer interface in all cases and, further, that the magnitude of 
these fluctuations is largely unaffected by flow direction (figure 5.10). 


5.4.3 Comparison of System Results at 33,5 rpm 

A comparison of the bulk velocity profiles at a rig speed of 33.5 rpm is 
shown in figure 5.13. The piston leakage effects observed in this case are 
different from those observed at 11.3 rpm (figure 5.2) since no net leakage is 
apparent; there may even be a net mass gain. Figures 5.16 and 5.17 reveal 
pressure variations about the ambient pressure of 0.98274 bars (as do figures 
5.6 and 5.7 at 11.3 rpm). Hence, a possible mechanism for the absence of net 
cyclic leakage at 33.5 rpm is that pressure differential driven leakage across 
the piston in one flow direction is canceled by mass gain when the flow 
reverses. However, this is not a totally satisfactory explanation in light of 
figures 5.6 and 5.7 (reversing pressure gradient with net leakage (figure 
5.3)). A more probable hypothesis is that at 33.5 rpm, the period during 
which leakage flow exists is significantly smaller than that at 11.3 rpm. It 
may also be true that at higher speeds, the piston seals function better as a 
result of increased seal temperatures caused by greater friction. The net 
conclusion, therefore, is that piston leakage is an important boundary 
condition for the test section simulation and it requires more sophisticated 
experimental and modelling attention. 

Ignoring the irregularity in the experimental bulk velocity profiles 
alluded to earlier (particularly at about 122°), the predicted and 
experimental bulk velocity profiles of figure 5.13 are generally in phase. 

The phase agreement of the one -dimensional profile is better than that of the 
two-dimensional profile, which may be ascribed to differences in dissipation 
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FIGURE 5.14 SYSTEM AXIAL VELOCITY PROFILE FROM 1-D CODE: 33.5 rpm 
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FIGURE 5.15 SYSTEM AXIAL VELOCITY PROFILE FROM 2-D CODE: 33.5rp 













CO 

CD 

o 

cu 

^T 

CO 

03 

CZ 3 


r^ 

CO 

in 

TT 

03 

03 

03 

03 

03 

03 

03 

O 

O 

O 

o 

O 

O 

O 


o 

CO 

o cn 
o 

o cn ^ 

. CO ’ri 

o cn o 

co 

cn .. •«-* 

x 

.. < 

XXX 

£ -2 
^ CD L 
e (o 
o O ▼-» -Q 

o* o — ' • 2 

OJ o Z CD UJ O) 

-o • o cn cc 

*— o IH I 3 o 

uj i— cn 

_j . . hh .. cn .. 

0 z cn z uj z 

z W O M CE M 

< X CL X Q- X 


C /3 

cn 

cn 

n 

i-i 

HH 

X 

X 

X 

< 

< 

< 

X 

>- 

Kl 


120 


x 


FIGURE 5.16 SYSTEM PRESSURE PROFILE FROM 1-D CODE: 33.5rpm 
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FIGURE 5.17 SYSTEM PRESSURE PROFILE FROM 2-0 CODE: 33.5rpm 


modelling. As noted above, the pseudo-steady-state friction factor ba ®®^ 
calculation method used in the one -dimensional code produces more dissip 
than the two-dimensional code. The reduction in simulated dissipation at 33.5 
rpm (peak Reynolds numbers greater than 131,000) caused by neglecting 
turbulence in the two-dimensional model is thus manifest on a sys em 
manifestation absent at 11.3 rpm (see figures 5.2 and 5.3). 

The system axial velocity profiles depicted in figures 5.14 and 5.15 are 
identical in shape: the reduction in axial discretisation is apparent in the 
profile generated by the two-dimensional code. Comparing thesefiguresw t 
those generated at 11.3 rpm (figures 5.4 and 5.5), it is evident that the 
system axial velocity profile topology is independent of rig speed, which 
principally affects the peak flow velocity. This is consistent for an p 
ended system operating at atmospheric conditions. 

The system pressure profiles produced by the one- and two-dimensional 
codes at 33 5 rpm shown in figures 5.16 and 5.17 are similar in shape to their 
rniintemarts at 11.3 rpm given by figures 5.6 and 5.7, respectively. Th 
features of these latter figures are accentuated at 33.5 rpm particularly n 
the two-dimensional case (figure 5.17). The impact of lower bulk mean fluid 
temperatures simulated by the two-dimensional code in the test section i 
readily apparent. Figure 5.20 reveals a radial ^“Petature difference between 
the wall and the centerline exceeding 2° compared with .3 at 11.3 rpm This 
£ pril facie evidence of the importance of explicitly including turbulence 
phenomena in the simulation as well as the garnering °f sufficient 
experimental heat transfer data to investigate the validity of the two 
dimensional code thermal predictions. The jaggedness of figure 5.17 compared 
with its one-dimensional counterpart in figure 5.16 is notable the 

iaeeedness is more pronounced than that occurring at 11.3 rpm (figure 5./). 

»ble“«d pressure perturbations are a direct result of the temperature 
Itriations shown in figures 5.18 and 5.19 (which again are significantly 
ereater than at 11.3 rpm as reflected by figures 5.8 and 5.9). This i 
largely ascribable to the increased turbulence level, which produces peak RMS 
velocity fluctuations of 5.6 m/s at 33.5 rpm (figure 5.21) compared with only 

1.6 m/s at 11.3 rpm. 

Ignoring the unfiltered turbulence effects, further comparison of 
figures 5 18 and 5.19 shows that the basic temperature profile of the r g 
33 5 rpm is mora complex and less isothermal than at 11^3 rpeu A comparison 
of figures 5.8 and 5.18 (one-dimensional predictions) shows that the flow is 
not nearly as isothermal over the test section at 33.5 rpm as it is at 11.3 
rpm Of particular note is the relatively large drop in temperature between 
theheat exchangers and the test section. Part of this temperature gradient 
occurs as a result of the flow acceleration/deceleration in the nozzles where 
the reversible thermal/mechanical energy transfer (vdp) is greatest (see 
figures 5.16 and 5.17). This emphasizes the importance of generat ng i accurate 
boundary condition data for any two-dimensional simulation of the test 
section It also suggests that the rig instrumentation should be upgraded to 
include as complete a definition of transient test section boundary conditions 
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FIGURE 5.18 SYSTEM TEMPERATURE PROFILE FROM 1-D CODE: 33.5rpm 
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FIGURE 5.19 SYSTEM TEMPERATURE PROFILE FROM 2 D CODE: 33.5rpm 
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FIGURE 5.20 SIMUL. RADIAL TEMPERATURE PROFILE AT 1304mm: 33.5rpm 





as possible. This is extremely important when attempting to apply the rig 
test data to actual Stirling hardware, which generally has very different 
heater and cooler tube boundary conditions than those of the rig. 

Furthermore, it is important to ascertain the influence of boundary conditions 
on the test data since, particularly at elevated rig speeds, boundary 
condition effects may overwhelm the physical similarity between the rig and 
Stirling machine heat exchangers. 

At 33 5 rpm, the radial temperature profile of figure 5.20 is comparable 
in shape to that at 11.3 rpm (figure 5.10). Other than the radial temperature 
difference noted above, careful inspection shows figure 5.10 (at 11.3 rpm) to 
have irregularities at the points of flow reversal (0° and 180 , figure 5.5) 
which are absent at 33.5 rpm. The reasons for this are not clear and require 
further investigation. 

Many of the features of the RMS velocity fluctuation profile measured at 
11.3 rpm (figure 5.12) are evident at 33.5 rpm in figure 5.21. The centerline 
fluctuations appear to be relatively smaller in magnitude at 33.5 rpm while 
the degree of irregularity toward the tube wall apparently is not much 
affected by the rig speed. The monotonic relationship between rig speed and 
peak RMS fluctuation is also evident; a threefold increase in rig speed 
produces a peak RMS velocity fluctuation ratio of 3.5. This may prove useful 
in devising better friction factor and heat transfer factor correlations for 
use in one -dimensional simulations. 


5 . 4.4 Comparison of Test Section Results at 11,3 rpm 

Cyclic perspectives on the simulated and measured radial profiles of 
mean axial velocity at 11.3 rpm and at the center of the test section (probe 
position of 1304 mm) are given in figures 5.22 and 5.23. The absence of 
negative velocities in figure 5.22 is a result of hot wire anemometry yielding 
data that have no flow direction discrimination. Clearly the simulated and 
experimental profiles have different shapes, with the experimental profile 
being more rounded. As noted for figures 5.2 and 5.3, the peak velocities 
simulated and measured are not dissimilar, their difference being related to 
piston leakage effects. Mechanically induced flow irregularities are also 
apparent in figure 5.22 (as discussed previously). Of particular note in 
figure 5.23 is the slight occurrence of positive and negative flows (the 
shape velocity profile) existing at the flow reversal points of 0 and 180“ 
(visible at 0° and 360°). Bearing in mind the absolute nature of the 
experimental data, the same effect is discernable in figure 5.22 at 360“ where 
the absolute value of the velocity is greater toward the tube wall than at the 
centerline, which is not usually the case. 

A better comparison of the shape of the radial velocity profiles is 
depicted i* figures 5.24 and 5.25 which are superimposed in figure 5.26. At 
210“ to 275°, the simulated profiles are much flatter than their experimental 
counterparts although their terminating values at the tube centerline and at 
about 18.5 mm are comparable. At 200“, both experimental and simulated 
profiles have a negative radial gradient in the core flow region although the 
simulated gradient negativity is more pronounced. 
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FIGURE 5.21 EXP. RADIAL PROFILE OF AXIAL RMS VEL. FLUCT. AT 1304mm: 33.5rpm 






FIGURE 5.22 EXP. RADIAL PROFILE OF MEAN AXIAL VEL. AT 1304mm: 11.3rpm 
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FIGURE 5,23 SIMUL. RADIAL PROFILE OF MEAN AXIAL VEL . AT 1304mm: 11.3rpm 
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EXPERIMENTAL RADIAL PROFILES OF MEAN AXIAL VELOCITY 
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Figure 5.25 



COMBINED RADIAL PROFILES OF MEAN AXIAL VELOCITY 

RIG SPEED: 11.3 rpm PROBE POSITION: 1304 mm 
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An alternate explanation for the disagreement between the shape of the 
simulated and measured profiles could be sought in a mismatch of experimental 
and simulated test section boundary conditions. In other words, if the 
experimental flow at the test section / nozzle interface contains an embedded 
radial gradient of axial velocity, then the roundness of the experimental 
profiles could be a distorted version of the embedded boundary gradient. The 
absence of an embedded boundary gradient in the simulation (slab flow 
conditions are assumed at the test section entrance) would account for the 
observed discrepancy. However, figures 5.29 and 5.30 do not substantiate this 
hypothesis since the velocity profiles at a probe position of 609 mm from the 
entrance are much flatter than those at the test section midpoint. This 
behavior is consistent with the existence of slab flow at the test section 
entrance under experimental conditions. Thus it appears necessary to seek an 
explanation for the observed discrepancies in terms of other factors such as 
the exclusion of turbulence from the simulation. 


5 . 4.5 Comparison of Test Section Results at 33.5 rpm 

The difference in shape noted for figures 5.22 and 5.23 at 11.3 rpm is 
more pronounced at 33.5 rpm as shown by figures 5.31 and 5.32. In this case, 
the roundness of the experimental radial profile of mean axial velocity 
(figure 5.31) is sharply contrasted against figure 5.32, whose flatness at 

33.5 rpm is more marked than that of figure 5.23 at 11.3 rpm, particularly at 
the boundary layer interface. The reversed boundary flow (relative to the 
direction of the flow core) at 0 and 360° is less than that revealed by 
figures 5.22 and 5.23; indeed, its existence in figure 5.31 is arguable. 

The shapes of the radial profiles at discrete crank angles are shown 
separately for the simulated and experimental data in figures 5.33 and 5.34 
respectively and together in figure 5.35. The reduction in roundness of the 
simulated profiles at 33.5 rpm compared with those at 11.3 rpm (figure 5.24) 
is evident. Figure 5.35 shows that the turbulence -induced radial flow effects 
are more apparent at 33.5 rpm than at 11.3 rpm. This is consistent with an 
increase in peak measured RMS fluctuation from 1.6 m/s at 11.3 rpm to 5.6 m/s 
at 33.5 rpm. The simulated peak radial flow velocity shown in figure 5.36 is 
still small at -0.011 m/s; two orders of magnitude smaller than the turbulence 
fluctuations . Thus momentum conservation dictates that including higher 
negative radial velocities in the simulation would coerce the concurrence of 
figures 5.33 and 5.34. In this regard, it may also be noted that the 
volumetric flow rates produced by integrating the experimental and simulated 
axial velocity profiles with respect to radius would be similar (allowing for 
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FIGURE 5.28 SIMUL. RADIAL PRESSURE PROFILE AT 1304mm: 11.3 rpm 





FIGURE 5.29 EXP. RADIAL PROFILE OF MEAN AXIAL VEL. AT 609mm: 11.3rpm 
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FIGURE 5.31 EXP. RADIAL PROFILE OF MEAN AXIAL VEL . AT 1304mm: 33.5rpm 
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FIGURE 5.32 SIMUL. RADIAL PROFILE OF MEAN AXIAL VEL. AT 1304mm: 33.5rpm 
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FIGURE 5.37 SIMUL. RADIAL PRESSURE PROFILE AT 1304mm: 33.5rpm 



piston leakage effects) as expected from the bulk velocity profiles shown in 
figure 5.13; in other words, the experimental and simulated results are 
consistent on a mass conservation basis. 

The absence of any discernable radial pressure gradient in figure 5.37 
is in accordance with the smallness of the radial velocities displayed by 
figure 5.36. The behavior revealed by figure 5.36 (as well as by figure 5.27 
at 11.3 rpm) does not appear to be in conflict with the experimental data, 
namely, a radial flow reversal corresponding with the axial flow reversal at 
0° and 180° (see figures 5.23 and 5.32). The extent to which this radial flow 
reversal is important in the laminar/turbulent flow transition process should 
be an important ingredient of future simulation and experimental 
investigations . 


5.4.6 Simulated Cyclic Performance 

An overview of the simulated cyclic performance of the test rig is given 
by table 5.3. The 'net kinetic energy conversion' is the sum of the 
irreversible kinetic/thermal energy conversion (or dissipation) and the 
isentropic heat generation rate (or v*VP 'work'). A positive value indicates 
that over the cycle, there is a net conversion of kinetic into thermal energy; 
a negative value indicates the reverse. Negative 'indicated work done' 
denotes that work is done on the fluid by the piston. 


Table 5.3 Base case cyclic energy balances. 


Engine 

Speed 

(rpm) 

Simul . 
Code 

Net Heat 
Input 

(j) 

Net Heat 
Output 
(J) 

Net Kinetic 
Energy Con- 
version^) 

Indicated 
Work Done 
(J) 

Energy 
Balance 
Error (%) 

11.3 

l-d 

3.3101 

23.1855 

-0.0165 

-19.8918 

0.000065 


2-d 

1.2399 

13.5442 

3.887 

-8.4173 

0.000241 

33.5 

l-d 

20.8935 

155.8568 

-0.0259 

-134.9901 

0.000534 


2-d 

6.5336 

45.4136 

-1.4571 

-40.332 

0.006644 


The magnitude of the difference in dissipation modelled by the one- and two- 
dimensional codes is manifested by the large difference in indicated work 
predictions. At 11.3 rpm, the one -dimensional code predicts an indicated work 
input 2.4 times greater than the two-dimensional code,, while at 33.5 rpm, the 
one -dimensional prediction is 3.3 times greater. The net heat output required 
to balance the indicated work input follows this trend while the net heat 
input is adjusted to close the energy balance. The two-dimensional net 
kinetic energy conversion values are large relative to their one -dimensional 
counterparts. This also is consistent with the under-prediction of 
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dissipation in the two-dimensional code owing to the neglect of turbulence in 
the test section. All the energy balance errors are less than .01% at worst 
with the worst error being achieved after 3 cycles from rest. 


5.4.7 Closure 

The salient conclusion to be drawn from the comparison of the baseline 
simulated and experimental results is that the two-dimensional simulation 
should include an adequate turbulence model in order to match the experimental 
data. The major physical impact of the measured turbulence appears to be a 
large increase in the radial advection fluxes in the test section compared 
with those generated under assumed laminar flow conditions. The effectiveness 
of the ' turbulence model* used in the one -dimensional code (that is, the 
pseudo- steady- state friction factor correlation) in replicating experimental 
conditions cannot be judged definitively due to the absence of adequate 
measurements. At the very least, such measurements should include ambient 
pressure and temperature as well as a transient pressure drop profile across 
the test section. 

However, the comparisons performed appear to confirm that the use of the 
pseudo -steady -state friction factor approach to dissipation modelling is at 
least an approximation that produces the correct qualitative trends in the 
simulated results. The importance of correct boundary condition modelling as 
an integral part of the two-dimensional test section simulation is also 
highlighted. Hence, it is doubtful whether generic two-dimensional codes that 
just simulate the test section can be expected to yield good quantitative 
answers which may be generalised to Stirling machine heat exchangers with 
arbitrarily complex boundary conditions. 

In terms of practical design application, these results point to the 
continued efficacy of a one-dimensional system analysis of Stirling machine 
fluid dynamics. The simulation/experimental result comparison shows that 
there is promise for modifying the pseudo -steady- state approach to reflect 
adequately the turbulence effects which, based on these data, seem to be a 
major cause of some of the anomalies in one -dimensional simulation so often 
cited in the literature. 


5.5 COMPARISON CASE SIMULATION APPLICATION 

The METR experimental parameters chosen for the development of the 
turbulence model are selected to represent the heater of the Space Power 
Demonstrator Engine (SPDE) operating at relatively high power. These 
parameters are shown in table 5.4. 
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Table 5.4 Comparison case experimental parameters. 


Operating point: MTI Experimental test no. 42 of 9/11/86 

Experimental indicated power: 17.9 kW 

Max. Reynolds number within the heater over a cycle: 13,985.3 

Maximum Valensi number within the heater over a cycle: 91.3 

Minimum Valensi number within the heater over a cycle: 71.7 

Maximum heater bulk Reynolds number over a cycle: 11,706.4 

Maximum heater bulk Valensi number over a cycle: 86.9 

Minimum heater bulk Valensi number over a cycle: 72.2 


where : 

^Re^bulk “ d Z(gV) / Z(nV) (5.1) 

(^Va^bulk “ 2(W) / Z(fiV) (5.2) 


Two sets of Reynolds and Valensi numbers are defined. One set 
corresponds to the maxima and minima occurring within the heater while the 
second set defines 'bulk' values, namely those obtained from volume averages 
over the entire heater (equations (5.1) and (5.2)). In terms of the METR 
operating philosophy, the rig parameters were chosen to reflect the bulk 
Reynolds and Valensi numbers with the latter taken as the mean of the minimum 
and maximum values. 

The adjusted rig parameters for the comparison case are noted in table 
5.5. The chief differences from the baseline case occur in the expansion 
space, conical transition piece, and flexible joint parameters. The test 
section length is reduced by 323 mm and a ' transition length' within the test 
section has been defined to contain the transition from a one- to a two- 
dimensional spatial discretisation. Previously, this transition was contained 
within the nozzle and hence occurred under axially varying flow area 
conditions. Placing the transition in the constant diameter test section is 
analytically preferable although no discernable numerical differences can be 
detected . 


5.6 THE TURBULENCE MODEL 

The only difference between the equation set used in the baseline and 
comparison cases is that the Reynolds stress tensor and turbulent heat flux 
terms have been reintroduced into the momentum and thermal energy equations 
respectively (equations (2.42.1) and (2.43.1)). 

The six restrictions constituting the turbulence model are listed in 
section 2.3. However, of principal concern in this discussion is the 
selection of a model fulfilling the stipulation of restriction V. Currently, 
the relevant turbulence terms in the momentum and thermal energy equations are 
represented by pseudo- steady- state friction factor and heat transfer 
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coefficient correlations. The purpose here, therefore, is to replace these 
correlations with a more fundamental turbulence model. 

The selection of the turbulence model is based upon the following 
criteria: 

a. General coordinate and dimensional invariance . 

Essentially, this means that the form of the model should be 
independent of the spatial discretisation scheme and system of 
units used. Casting the model in tensor form and using 
dimensionless constants meets this criterion. 

b. Full realizability (Lu78). 

In particular, this requires that all the components of the 
turbulent kinetic energy are non-negative for any possible flow 
condition. Further, full realizability means that the turbulence 
model must be capable of evolving positive turbulent kinetic 
energy from a rest state (that is, zero momentum) without the 
necessity of defining an artificial rest-state positive turbulence 
field. This is particularly important under oscillating flow 
conditions during which repeated laminar- to- turbulent transitions 
occur, since under these circumstances, numerically ill- 
conditioned models show steady- state solutions that are dependent 
on the initial conditions. However, physical realizability alone 
is insufficient in a numerical context since models that are 
physically realizable can still produce negative turbulent kinetic 
energies due to the errors inherent in a discrete numerical 
process, particularly in laminar- to -turbulent transitions. Thus 
the model must also be numerically realizable in the sense that 
its numerical or discrete analog will also yield non-negative 
turbulence kinetic energies under all flow conditions. 

c. Material frame indifference . 

This requires that the form of the Reynolds stress term (equation 
(2.42.2)) in a two-dimensional turbulence field be the same in 
all non- inertial frames of reference which can undergo arbitrary 
time -dependent rotations and translations relative to an inertial 
framing' (Sp87). Basically, this implies that if the Reynolds 
stress is dependent upon the mean velocity gradient tensor Vv, 
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Table 5.5 METR comparison simulation input parameter set. 


Expansion Space 


Piston stroke(mm) 

: 251.46 

Piston end cap thickness (mm) 

: 0. 

Piston end cap taper(deg) 

: 90. 

Clearance length at piston TDC(mm) 

: 67.9 

Cylinder diameter(ram) 

: 127. 

Angular velocitv(rpm) 

: 33.75 

Conical Transition Piece 


Axial length(mm) 

: 219.1 

Cylinder mating diameter(mm) 

: 127.0 

Flexible joint mating diameter(mm) 

: 127.0 

# of control volumes 

: 3 

Flexible Joint 


Cylinder side axial length (mm) 

15.875 

Cylinder side diameter(mm) 

152.4 

Cylinder side # of control volumes 

2 

Nozzle side axial length(mm) 

73.025 

Nozzle side diameter(mm) 

139.7 

Nozzle side # of control volumes 

2 

Spacer 


Axial length(mm) 

63.5 

Cylinder side flange length(mm) 

3.175 

Nozzle side flange length(mm) 

6.35 

Diameter (mm) 

139.7 

Cylinder side flange diameter(mm) 

127.0 

Nozzle side flange diameter (mm) 

127.0 

# of control volumes 

3 

Heat Exchanger 


Axial length(mm) 

50.8 

Cylinder side flange length(mm) 

9.525 

Nozzle side flange length(mm) 

12.7 

Casing diameter(mm) 

127.0 

Hydraulic diameter (mm) 

2.1824 

Flow area ratio 

.556 

# of control volumes 

4 

Nozzle 


Size (1-small, 2-large) : 

1 

# of half-nozzle control vols 

3 

Test Section 


External diameter(mm) : 

50.8 

Axial length(mm) : 

2286.0 

Transition length(diameters) ; 

2 

# of axial control volumes (4 minimum) : 

11 (2-d) 

# of radial control volumes (2 minimum) 

8 

Wall conductivity(W/ra-K) 

1.0211 

Wall density(kg/m A 3) : 

2466.84 

Wall specific heat caoacitvf J/ke-IO 

753.624 

Working Fluid 


Working fluid (1-Air. 2-Helium. S-Hvdroeen’) : 

1 
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then Vv must always be non- zero In a two-dimensional turbulence 
field and a transformation tensor K must exist such that: 

V*v* - Vv + dual K 

where * indicates the target non- inertial frame of reference. 

A fourth (although non-rigorous) selection criterion is that the model 
should be amenable to the induction of and the response to temporal 
discontinuities. This arises because of the nature of the turbulence 
triggering occurring in oscillating flows. An example of such a discontinuity 
is evident in figure 5.38 (Se88) , which is the closest test point in Seume's 
data set to the conditions of table 5.4. The abrupt onset of turbulence at 
255° and to a lesser extent at about 20° as well as the decay of turbulence at 
303° point to discontinuous phenomena. In this case, Seume postulates that 
the turbulence is triggered advectively and decays rapidly after passage of a 
turbulent fluid slug. In contrast, Hino et al. (HK83) concluded from their 
oscillating flow experimental data taken in a rectangular duct that turbulence 
is generated 'explosively' at the onset of flow deceleration, which suggests 
that local fluid acceleration is the trigger mechanism. Hence in this case as 
well, the turbulence triggering is a discontinuous event. 

The approach decided upon for implementing these empirically observed 
discontinuities in the model is via a foldback function formulation. This 
formulation has the advantage of introducing and controlling discontinuities 
by continuously changing selected parameters, an ideal characteristic for a 
numerical implementation. An example of a simple foldback function is given 
in figures 5.39 to 5.42. Consider a parabolic function of the form: 


y - ax 2 + bx - f(x) 
Hence, initially if: 
yi — aXQ 2 + bxQ — f (x 0 ) 


then, if yi is folded back into f(x), that is: 


72 “ f(7i) “ a 7i 2 + b Yx 


(5.4.2) 


equations (5.4) represent a foldback function. This function is depicted in 
figure 5.39 where the dashed line represents the foldback process as the 
output of the function is reflected by the 45° line back into the function 
input. Hence if a - -1 and 6 - 3 , the function is 'ordered' since the 
solution converges to a unique value. Setting b - 4 produces the effect of 
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figure 5.40, where the solution converges to a 'stable' condition by tracing a 
closed coordinate path. Increasing b to 5 yields a 'transition' condition in 
which the solution gradually degenerates into chaos (figure 5.41). Finally, 
figure 5.42 shows a fully chaotic condition when a - -.5 . Hence a continuous 
change in a and/or b produces a discontinuous change in the stability of the 
foldback function so that the chaos 'evolves' in a fashion simulating that 
observed empirically without creating numerical anomalies. 

The foldback function methodology invoked in this analysis was selected 
initially because it is ideally suited to the iterative structure of the 
numerical algorithm used (section 2.5), which strongly motivated its initial 
selection. Yet it must be emphasized that other classes of 'chaos' models 
(such as fractal methods) are potentially useful if not superior for modelling 
laminar -to -turbulent (and turbulent- to -laminar) transitions. Time 
constraints, however, did not permit the investigation of alternates to the 
foldback function methodology chosen. 


5.6.1 Selection of the Turbulence Model 


Consider the Reynolds or turbulent shear stress term in the momentum 
equation (2.42.1): 


[V ]T (t) - - tv]g' [V] v ' 


In terms of the rigorous discrete volume analysis used as the basis of the 
simulation (Go87) together with the six restrictions defining the turbulence 
model the ideal choice of turbulence model would be a functional form 
yielding the correlation on the right-hand side (RHS) of equation (5.5) 
directly By a manipulation of the Navier-Stokes equations (Hi75), the exact 
differential transport equations for the Reynolds stresses may be expressed as 

(HL72) : 


D(g£v'.)/Dt - -(vjgi- d Vi /dx k + gjvjf 3vj/dx k) 

generation 


-2/i(3vj/3x k )(3vy3x fc ) + P'idv^/dxj + dv'j/dXi) 

dissipation redistribution 


-d[g' i v' j vj c - |i{3(v}vy/3x fc + P'(«jfcv< + S ik vp]/dx k 

diffusion 


where the subscripts i, j, and k imply the usual tensor summations. In its 
exact form, equation (5.6) is not soluble since the RHS contains a number of 
new correlations (for example, in the dissipation and diffusion terms) for 
which closed form determinations must be found, that is, the ubiquitous 
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closure problem. Typically, these closures are achieved (for example, see 
HL72, LR75, GE76) by specifying the unknown correlations in terms of algebraic 
functions of known correlations, the turbulence kinetic energy (fc) , and the 
rate of dissipation of turbulence kinetic energy (e), which are defined as 
follows : 


k - 0 . 5(v|) 2 


(5.7.1) 


e- t/(3v| /dxj) 2 


(5.7.2) 


Employing these techniques expands equation (5.6) into a set of six equations 
for the two-dimensional field simulated in the METR test section, namely, four 
equations for the Reynolds stresses and two for k and e , 
respectively. 

However, it is also possible to determine the Reynolds stresses directly 
from a k-e equation set using Boussinesq's eddy-viscosity concept, which 
assumes that the turbulent shear stresses are proportional to the mean- 
velocity gradients (analogous to the viscous stresses in laminar flows) or: 


T ij “ M (t) (dvj/3xj + 3vy3x i ) - 2kS i j/3 


(5.8) 


The turbulent viscosity p (t) can be obtained from the Kolmogorov- Prandtl 
expression (Ro84) : ' 


/i (t> - pK k 2 /e 


(5.9) 


Equation (5.8), however, suffers from a major disadvantage because it yields 
isotropic principal Reynolds stresses, that is: 


11 T JJ 



(5.10) 


This has the consequence of erroneously predicting uni-directional flows in 
non-circular ducts as well as producing substantial inaccuracies in the 
calculation of separated turbulent flows (such as that occurring over a 
backward facing step). Speziale (Sp87) has proposed solving this problem by 
replacing the linear form of equation (5.8) with a non-linear, quadratic form. 
He has solved several test problems successfully using this approach. 

Therefore, in the context of the METR analysis, since a closure model of 
the two -parameter or k-e type is in any event necessary for implementing the 
desirable full Reynolds stress model, in terms of the imposed temporal and 


155 



fiscal constraints, it was judged prudent to limit the turbulence model in the 
METR simulation to such a two-parameter model, at least as a first step. 
Equations (5.8) and (5.9) are thus used to generate the necessary Reynolds 
stresses. Clearly, the limitations of equation (5.8) must be acknowledged, 
but, fortunately, they are of second order significance for the constant area, 
axi- symmetric METR test section to which the turbulence model is being applied 
in two dimensions. Hopefully, the experience gained with the k-t class of 
two-parameter turbulence model as applied to oscillating flows will provide a 
basis for future work leading to the implementation of a full Reynolds stress 

model. 


5.6.2 The k-£ Turbulence Model 

Examination of the classic k-t turbulent model provides a convenient 
starting point for the model selection process. A general form of the k-t 
model in common usage is referred to as the 'low Reynolds number version 
which is purportedly valid throughout the laminar, semi -laminar, and turbulent 
regions of the flow (LS74) . Usually, these characteristics are attributed to 
the model under steady- state flow conditions. It is not clear from the 
literature whether the model is intrinsically capable of handling general 
laminar -to -turbulent transitions (including those where the fluid is 
accelerated from rest), although local turbulent -to -laminar transitions have 
been modelled (JL72). The general form of the k-t model also includes wa 
correction terms not only to account for near-wall effects but also to allow 
setting £ - 0 at the wall (JL72). In terms of the turbulence model 
restrictions, a tensor notation version of the low Reynolds number k-t model 
suitable for compressible fluids has been adapted from LS74 and is given by 
the following differential transport equations: 


pDk/Dt - V.{(/i+/C 1 /i <t) )Vk) + 2/i (t) (D:Vv) - 2 MVk 1 ' 2 ) 2 - P< 

diffusion production wall correction dissipation 


(5.11.1) 


pDe/Dt - V.{(p+Jf 2 P (t) )^) + *3 e2p (t) (°: ^v)/k 

diffusion production 


(5.11.2) 


- K^pt 2 /k + 2pp (t) (^2* 7v ) 2 

destruction wall correction 


These equations share a common form having diffusion, production, 
destruction (or dissipation), and wall correction terms on the RHS. Of 
particular note is the form of the production terms which are defined by the 
scalar protect of the deformation tensor (equation 5.11.4) and the velocity 
gradient. As this product is non-negative under all flow conditions the 
creation or evolution of k and £ is dependent on the value of >, itself a 
function of k and £ (equation 5.11.7). Hence, the system of equations is 
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indeterminate even though it is closed. The balance of the equations filling 
out the complete model implementation are given by: 


T (t) - - 2pkl/3 

D - 0.5{Vv + (Vv) T ) 

-Si V j - T?] 

gp" “ H^HdT/dx^/Np' 

/* (t) - Ksptf/e 

K i ~ 1 K z - 1/1.3 k 3 - 1.44 

g* “ 1.92(1 - 0.3exp(-tfj^ 2 ) } 

K 5 - 0.09exp{ -2. 5/(1 + tf$2/50) 


(5.11.3) 

(5.11.4) 

(5.11.5) 

(5.11.6) 

(5.11.7) 

(5.11.8) 


N He ~ Pk 2 /pe 


(5.11.9) 


Equation (5.11 6) defines the method by which turbulent enthalpy transport 
(equation ( 2 35) ) is determined from the model while the constants K< and K 5 
ependent on the turbulent Reynolds number (equation (5.11.9)). 

_ j i Examination of equations (5.11.1) and (5.11.2) reveals that the k-e 
model satisfies the first selection criterion (general coordinate and 
dimensional invariance) while equation (5.11.3) satisfies the third criterion 
(material frame indifference). However, the second criterion (full 
realizability) is not satisfied. V 

F irst, equations (5.11.7) and (5.11.9) are indeterminate if e - 0 and 

equation (5.11.2) is indeterminate if k - 0. Hence the equations are 

n rinsically incapable of evolving a positive turbulence while either k or e 

nIp Z ? r ?M S* us ^ al methodol °gy for dealing with this situation is simply to 

P i« ni K 4 ?^ iSe the k Snd C fields to values that are compatible with the^^ 

^1 ty 1 H r K iteri °? : namely * the sourc e terms in equations (5.11.1) and 

(5.11.2) should be positive, or: ' 


(2/i (t) (D:Vv) - 2p(Vk 1 ' 2 ) 2 


P^init > 0 


(5.12.1) 



(5.12.2) 


UC 3 e2M (t) (D:Vv) - K l ,pe 2 /k + 2 M^ (t) (V I -Vv) 2 } init > 0 


In this sense, the k-e model may be termed 'weakly' realizable. 

However, in oscillating flows where k and « can be zero in the flow core 
as well as in the laminar boundary layers, such initialisation techniques 
mitigate against the prediction of strict laminar- to -turbulent transitions 
(that is, from k - 0 to k > 0) by the model per se. Furthermore, from a 
numerical perspective, there does not appear to be a physically credible way 
of ensuring that the k source terms (LHS of equations (5.12 1)) will always 
interact with the transport terms in such a way as to unconditionally produce 
positive values of k when considering the model from a full realizability 
perspective . 

As a result of its system indeterminacy, the k-e model is suitable for 
implementation as a foldback function, so fulfilling the fourth selection 
criterion. This may be readily achieved by casting equations (5.11.1) and 
(5 11 2) in implicit or advanced time form, thus rendering them nature y 
amenable to iterative solution. There are then a number of possibilities 
(notably in the generation terms) for introducing foldback parameters for 
continuously controlling the triggering and evolution of the turbulence. This 
capability opened the possibility that full realizability could be achl ®^ 
for the model when cast in foldback function form, thus warranting a numerical 

investigation. 

After much tedious experimentation, it eventually became clear that the 
realizability problem proved to be the main detractor in obtaining anything 
resembling the type of laminar- to -turbulent transition observed empirically y 
Seume as fhown in figure 5.38. The model could be made to 'work , that is, 
portray discontinuous transitions, by treating the flow as being turbulent ^ 
throughout the cycle (that is, k and e > 0 for all t) simply by setting p 
throug y . . (taken to be 2,000). Further difficulties were 

encounter* wifSfS'TK'&SSi* grid framework of the discrete voluce analysis 
(eight control volumes radially) as a result of the wall correction terms in 
equations (5.11.1) and (5.11.2). The spatial discretisation near the wall 
strongly determines the influence of these terms, thus violating the 
scalability inherent in the discrete volume analysis. Alternate wall 
correction procedures such as the logarithmic wall function method (LS74) were 
.Iso tested P with mixed results. Such empirical procedures by definition are 
not spatially scalable or generalized. 

These experiences prompted a search of the literature for alternate two- 
parameter turbulence models which could meet the full realizability criterion 
and also avoid the necessity of special wall functions. Accurate prediction 
of steady- tate turbulence or turbulent - to - laminar transition phenomena alone 
(commonly used for validation) was not regarded as a qualifying criterion. 
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5.6.3 The Wilcox k^u Turbulence model 


Wilcox (Wi88) has proposed a k-u two -parameter turbulence model which, 
when implemented in a foldback function, nominally comes much closer than the 
k- e model to achieving full realizability (which appears to be necessary for 
modelling discontinuous transitions). Although still not possessing complete 
system determinacy, the Wilcox k-u model does appear to offer a better 
opportunity of meeting the full realizability criterion than most of the two- 
parameter models reported (Ro84, TL72) . It may be noted that Spalding also 
published a report on the k-u model (Sp72), but it is not known to what extent 
the Spalding k-u model differs from the Wilcox k-u model, u is referred to as 
the specific dissipation rate and Is related to e via: 


w - e/Kk 


(5.13) 


After a series of numerical applications, Wilcox claims that the k-u 
model is 'much more accurate than comparable two-equation models' although 
such a claim is, of course, disputable. The Wilcox k-u model does not include 
any wall correction terms since the model can be integrated through the 
viscous sub-layer at the boundary. Hence these characteristics make the k-u 
model an ideal candidate for implementation in a discrete volume analysis. 

terms of the turbulence model restrictions, a tensor notation version 
of the k-u model given in Wi88 may be expressed by the following differential 
transport equations: 


pDk/Dt - V.{(/i+^2M (t> )V/c} 

diffusion 


+ (T (t) : Vv) 

production 


- K 1 pu>k 

dissipation 


(5.14.1) 


pDw/Dfc - V* { (/i+K 3 fi (t) )Vu) + K 5 u(T (t) :Vv)/k - K 6 pu 2 

diffusion production destruction 


(5.14.2) 


These equations are similar in form to their counterparts in the k-e model 
(equations (5.11.1) and (5.11.2)) except that the wall correction terms are 
absent. Another notable difference is that the production terms are expressed 
in terms of the Reynolds stress tensor rather than the deformation tensor. By 
an argument similar to that invoked previously, the k-u equation system is 
also indeterminate in explicit form. The balance of the equations and 
constants comprising the complete model are given by; 


T <t) _ 2/x (fc) {D - (Vv) 1/3 ) - 2pkI/3 


(5.14.3) 


D - 0 . 5 { Vv + (Vv) T ) 


(5.14.4) 
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.giu • * f CM 
8i v j T ij 


(5.14.5) 

- K^k/u 


(5.14.6) 

g ' i T > - j* <t, (ar/ax 1 )/NjJg 


(5.14.7) 

JCi — 0.09 K z - 0.5 

* 3 - 0.5; 

(5.14.8) 



JC A - 1 K 5 - 5/9 

* 6 “ 3/40 



In this case, none of the constants show a turbulent Reynolds number 
dependence. However, the methodology employed by Wilcox to derive the model 
constants (equation 5.14.8) is open to question since it relies partly on 
numerical evaluation (as opposed to being entirely based on experimental 
data) . This has particular significance in the light of the scalability 
requirements of a discrete volume application. 

The foldback function formulation of the k-u> model is achieved via the 
following manipulation. Substituting equation (5.14.6) into equation (5.14.3) 
yields : 


T (t > - 2K k pk(V - (Vv)I/3}/w 2picI/3 


(5.15) 


Substituting equation (5.15) into the production term of equation (5.14.2) 
produces : 


pIXo/D t - V.{(^ 3 /i (t, )Vc - * 6 pa> 2 + 2 pK 5 [ [K*{D- (Vv)I/3) - wI/3] :Vv ] (5.16) 


Casting equation (5.16) into advanced time (or implicit) form by discretising 
the LHS and rearranging the source (production and destruction) terms yields . 


pw[ 1/At + t (22C 5 I/3) :Vv) + K 6 w] - V* { (p+^ 3 p {t) )Vw) 

- 7« S /At + 2pK^K 5 [ { D - (Vv)I/3) :Vv] (5.17) 

Equation (5.17) is fully realizable. It is also determinate if in 
equation (5.14.6): 
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w — 0 


A* (t) - 0 


(5.18) 


Notice then that as long as the velocity gradient Vv is non-zero, equation 
(5.18) will evolve positive « even when w - 0 regardless of the value of k. 
This implies a much different physical description of turbulence than the k-e 
model, which will only evolve positive e when k > 0. A determination of which 
model is physically 'correct' under oscillating flow conditions (particularly 
during laminar- to- turbulent transitions) requires a level of analysis beyond 
the scope of this project and a very detailed experimental oscillating flow 
data base, which is not yet available (at least from the METR) . Nevertheless, 
Wilcox has demonstrated that under steady-state conditions and during 
turbulent -to -laminar transitions, the k-u model is at least as physically 
'correct' as the k-e model. Hence for the sake of expediency, it will be 
assumed here that the physical implications of equation (5.17) are valid on 
the basis that the production term of equation (5.17) is a logical consequence 
of equations (5.14). 

Casting equation (5.14.1) in implicit form yields: 


pk(l/At + - V*((p+* 2/1 <t>)V/c} - pfcy At + (T (t) :Vv) 


Equation (5-19) is also realizable (not as fully as equation (5.17) though) 
but it still suffers from the system indeterminacy of the k-e model in that 
positive k can only be evolved if k > 0. Hence some additional empiricism is 
necessary to trigger the creation of turbulence. 


5 - 6 - 4 ft-oj Model Enhancem ent and Foldback Function Parameters 

Dropping the implicit notation for the sake of clarity and transforming 

equations (5.17) and (5.19) into generalised discrete volume form (see chapter 
2) yields: 


d ([V]A: M ( s))/dt - _ tv n )*<( (fcVn ,g- ttV p v n(£r) ) . - n )dA 

J A ms) 

+ *3 [V]/ i<t) ) (V [Vn j]c* -n)dA (5.20) 

A ms ) 


^(S)(* 2 (tv]T (t) : V [tv] v) - 


^1 [tvjp [V] w [V]^} 


and: 
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_ [V n ] w ( ([tv n ]g-ttv n )P V n(S>) * 
J ^n(S) 


_ ([tv n ]M + *4 m ^ (t) )(^cv n i«*-n)dA 

J^n(S) 


(5.21) 


+ ^ 6 M (s) [2{ ItV] D - (V. [tV jV- m w)I/3):V ltvJ v] 


- *7 W (S) [V) w2 


for the k and « transport equations, respectively. Note that the constants K 
in these integral equations do not correspond exactly to those used in the 
equations of section 5.6.2. 

The remaining equations and constants constituting the implemented 
version of the k-u model will be listed first and discussed thereafter. Two 
turbulence triggering modes are investigated and are distinguished by 
referring to them as the 'discontinuous' and 'continuous' modes. The lirst 
block of equations and constants are common to both triggering modes and are 

given by: 


“ *5*10 [tV]P [V]* / tV] w 


(5.22.1) 

JC 5 - {1 - exp ( -r + / 26. 5) } 2 


(5.22.2) 

r + “ ([tv]P [tv] r wall)° 5 ( r wall * 

1 ^ 
| 

(5.22.3) 

K x - 0.09 K 3 - 1 *4 “ 0.5 

X 6 - (3 + 5K s )/9 Ky - 3/40 


(5.22.4) 

N $e ~ [tv ] P [V]V*i [tv]M [V) w 


(5.22.5) 

[tv) D - 0.5 {V (tv] v + (V [tV] v) T ) 


(5.22.6) 

[V] T (t) ~ 2 tV) /i <t) { [ t vi D ‘ (V[tv]V)I/3) 

- 2 ( [ -tvj P [V]*)I/3 

(5.22.7) 
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(5.22.8) 


[V]$ (t) ~ *1 [tV]P [V]^ tV] w 

[V n ]q (t) - Cp m^ (t) V [t v n] r/4y (5.22.9) 

The second block of equations and constants are differentiated by triggering 
mode as follows: 

discontinuous: K B - 5 [ log(l+K g ) A { l+exp(l-K g ) ) ]/3 

, (5.23.1) 

continuous: K B = min ( [ 5 [ log(l+K 9 ) A { l+exp(l-*g) ) ]/3 , l) 

discontinuous: K 9 - / (^critical 

continuous: * 9 - (Nj&-1)/(N^ Critical (5.2J.2) 


discontinuous 

continuous 



critical “ ^6 

critical " 25 


(5.23.3) 

*2-1 


00 

tN 


(5.23.4) 

*10 “ 1 

Re ) — transition 

*io - 1 


(5.23.5) 

- 0 

Re ) ^ (^R e ) transition 





Equations (5.22) and (5.23) embrace the foldback function parameters and 
include enhancements to the Wilcox k-u> model as discussed below. 

a. Near-wall correction methodology 

As alluded to previously, the Wilcox k-u model includes no near- 
wall correction terms. However, in contrast to the coarse radial 
discretisation typifying the discrete volume analysis used, a much finer 
spatial discretisation is necessary near the wall to achieve adequate 
physical accuracy with the basic model. This requirement substantially 
^■ nc t63.ses the computational cost of implementing a discrete volume 
analysis. Inclusion of an appropriate near-wall correction term, 
however, enables a much lower computational cost to be achieved while 
retaining scalability. Nagano and Hishida (NH87) investigated a range 
of existing near-wall correction terms as applied to a k-e turbulence 
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model and proposed an alternative which, they concluded, yields better 
performance than the existing methods. 

Using the k-w model, some investigations with both the Jones- 
Launder (JL72) and Nagano -Hishi da near-wall correction methods were 
performed with the latter method yielding better agreement with the test 
data at the chosen test point (section 5.5). Hence on this basis, the 
Nagano -Hishida method was selected pending the availability of further 
experimental data. The model consists of equations (5.22.2) and 
(5.22.3) which serve to reduce exponentially the magnitude of the 
turbulent viscosity as the wall is approached. A dimensionless distance 
(zero at the wall) defined in terms of the wall shear stress is used as 
the independent parameter. This ensures that the method is independent 
of spatial discretisation, so retaining the scalability characteristic 
of the discrete volume analysis. 

b. Constant modifications 

Other than the introduction of K$ (wall correction term) , K 2 
(continuous turbulence trigger parameter) and K 6 (foldback function 
parameter), with one exception, the constants are the same as those in 
the Wilcox k-w model. The exception concerns the turbulent diffusion 
term constant K 2 in equation (5.14.1) which Wilcox assigned a value of 
0.5 (equation (5.14.8)) based on numerical experiment. When applied to 
an oscillating flow, changing this constant to a value of unity (* 3 in 
equations (5.20) and (5.22.4)) yielded better agreement with the 
experimental data at the test point. Jones and Launder (JL72) and 
subsequent workers have also used a value of unity for the constant in 
the k equation turbulent diffusion term (equation (5.11.1)). 

c. Turbulence triggering 

As noted in the discussion of equation (5.19), the k-w model can 
only evolve positive turbulent kinetic energy if k is made positive 
independently. In general, it seems that two approaches to introducing 
this independence can be considered. 

As a first option, the advective triggering mechanism postulated 
by Seume (Se88) may be used to inject positive k into the two- 
dimensional test section of the METR. This is adequate for the general 
case provided that such positive k can be generated outside the test 
section, that is, in the one -dimensional duct between the test section 
and the heat exchangers. Once a positive fc thus has been established, 
the turbulence evolution may be controlled via a foldback function 
parameter such as the turbulent Reynolds number (equation ( 5 . 22 . 5 )). 

A second, more fundamental, option is based on the local 
acceleration triggering postulate of Hino et al. (HK83) . Such a 
postulate devolves to a triggering model based upon microscopic level 
entity (molecules or 'clumps' of molecules) dynamics in which random 
cross-stream or turbulence -inducing momentum arises when the microscopic 
viscous forces are exceeded by the inertial forces . This mechanism is 
depicted on a macroscopic level in contained flows by the usual Reynolds 
number transition criterion. Proceeding further, a microscopic entity 
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approach also would allow the discontinuous increase in turbulence 
observed by Seume (apparently caused by the passage of a turbulent fluid 
slug) to be naturally modelled so harmoniously blending Hino et al . f s 
and Seume 's observations. 

The microscopic turbulence triggering mechanism is considered 
preferable by the author; however, the project time and budget 
constraints mitigated against investigation of this concept. Hence, by 
default, it was decided to proceed with the advective triggering 
methodology even in the presence of concerns about its physical 
credibility. 

The advective triggering methodology is implemented by creating a 
pseudo- transport equation for k which is used only in the one- 
dimensional sections of the METR between the test section and the heat 
exchangers. All turbulence is assumed to be extinguished in the 
serpentine passages of the heat exchangers. The one -dimensional pseudo - 
transport equation for k is given by: 


d([ V ]/c Af (S ))/dt 


[V n ]fc { ([tV„] g 

J A ms) 


[ tv n 1 P v n(S))*' n ) dA 


(5.24.1) 


i f : 

[V]^ < [tv] T wall / [tv]P and N Re > 2000 

(5.24.2) 

then : _ _ 

[V]^ “ [t V) r wall / ctvj P 


Equation (5.24.1) just reflects advective transport of k while equation 
(5.24.2) initialises k to a floor value of the square of the friction 
velocity if the Reynolds number is greater than a transition value of 
2000. If k is greater than the floor value (irrespective of the 
Reynolds number) it retains the value computed using equation (5.24.1). 
This approach yields a fully realizable turbulence triggering method 
including laminar- to- turbulent transitions from the rest state without 
the necessity of any k or w pre-initialization. 

d. Foldback function implementation 

Control of the foldback function is achieved through constants K 2 
and K 6 (equations (5.23.4) and (5.22.4), respectively) together with 
their threaded substituent constants. 

A stability analysis of the k-u model showed that the w production 
term in equation (5.21) is an appropriate vehicle for controlling the 
stability of the foldback function. In particular, the stability 
boundary hypersurface (Go87) for the model can be intersected for a 

of Kq between 7/9 and 8/9 (all other constants remaining at their 


* 


I* 
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listed values). If Jt 6 < 7/9, the model exhibits chaos increasing in 
magnitude until a fully chaotic condition (figure 5.42) is reached for 
JC 6 = 3/9. For K 6 > 8/9, the model exhibits 'ordered' behavior (figure 
5.39). Thus, continuous variation of K 6 between 3/9 and 8/9 enables the 
k-w model discontinuous ly to evolve turbulence in a random fashion. 


The continuous variation of K e is achieved via K a (equation 
(5.23.1)), which varies as an inverse exponential between 0 and 1, 
asymptoting to a straight line of small positive slope (0.084) for Kg ^ 
1. JC 9 is defined as the ratio between the actual turbulent Reynolds 
number and a 'critical' value determined numerically on the basis of 
matching the available experimental data. 

e. Turbulence triggering modes 

The discontinuous turbulence triggering mode is implemented using 
the classic Reynolds number transition criterion for pipe flow as shown 
by equation (5.23.5). In essence, this zeroes /i (t> f° r Reynolds 

numbers less than the transition number, which is set at 2000. Note 
that K z is maintained at unity, thus placing no additional restrictions 
on the production of k while ) C ritical is assi g ned a relatively 

high value of 486, again, based on purely numerical considerations. 

The continuous turbulence triggering mode is achieved by allowing 
/i (t) to exist regardless of the value of k as long as w > 0 (equation 
(5.23.5)). However, the foldback function control parameter is modified 
slightly from its discontinuous form (equation (5.23)) and applied not 
only to the w production term in equation (5.21) but also to the k 
production term in equation (5.20). Production of k and w is thus 
independent of the bulk flow Reynolds number (K B has a floor value of 
unity in its continuous form). (Wj?e) critical is assi S ned a value based 
on the criterion that k is equal to the square of the friction velocity 
under 1 ami nar flow conditions. Physically, this implies a 'continuous' 
model of fluid turbulence in which all flows contain some level of 
turbulence except that under laminar conditions, such turbulence is 
small, stable or 'ordered' and does not increase in the presence of 
imposed stimuli. Atmospheric flows, for example, are often considered 
in this light. 

Thus, in summary, the k-w model is implemented as defined by equations 
(5.20) to (5.24) at the chosen experimental test point. In view of the 
foldback function methodology employed, the empirical nature of its 
constituent parameters, the uncertainty as to which triggering mode is 
physically accurate and the basic reservations expressed about the advection 
triggering methodology, this version of the k-w model should be viewed as 
entirely experimental and subject to even drastic revision. 


5.7 COMPARISON CASE EXPERIMENTAL DATA 

The experimental data made available for the comparison case simulation 
application consists of a half cycle of axial ensemble -averaged mean 
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velocities and root-mean-square (RMS) velocity fluctuations at 12 locations on 
one radial plane. The flow during the half cycle is in the negative 
direction, that is, from the open end toward the cylinder. The parameters 
defining the experimental data are listed in table 5.6. 


Table 5.6 Parameters for the comparison case experimental data. 


Engine speed: 33.75 rpm 

Probe position: 1676.4 mm from the cylinder end 

Dry bulb temperature: 25.5°C 

Wet bulb temperature: 15.8°C 

Ambient pressure: 981 bar 

No. of cycles for ensemble averaging: 500 

No. of probe radial locations: 12 


Since 500 cycles are used for computing the ensemble averages, in this 
respect, the data are felt to be adequate for validation purposes. The probe 
position is located toward the open end of the METR, so enhancing the yield of 
the comparison and validation process because any experimental flow 
asymmetries would stress the simulation to a greater extent. Unfortunately, 
since only negative half-cycle data are available, this aspect of the 
comparison cannot be accomplished. 

Unlike the base case data, no bulk mean velocities were provided for 
simulation calibration purposes. Hence these velocities are computed from the 
mean velocities provided using an area-weighted aggregation procedure. The 
resulting bulk velocity comparisons for the discontinuous and continuous 
triggering modes (henceforth referred to as the DTM and CTM simulations, 
respectively) are shown in figures 5.43 and 5.44. The simulated velocity 
profiles are identical showing that the triggering mode has no influence on 
the bulk velocities, as expected from a conservative set of transport 
equations. The experimental data reveal a small phase shift relative to the 
simulated data which is apparent at 0° and 180°. However, of major concern are 
the large amplitude, higher order harmonic fluctuations in the experimental 
data between 60° and 120°. Apparently, these fluctuations are related to 
piston/cylinder stiction and/or sealing effects. Nevertheless, their 
magnitude and shape has a major impact on the behavior of the flow turbulence, 
so complicating the comparison of the experimental and simulated k between 60° 
and 120° because the simulated k is not subject to the same production stimuli 
as the experimental k . Clearly, these fluctuations in the mean experimental 
data should be removed before any further comparison or validation work is 
undertaken . 


5.8 COMPARISON CASE RESULTS 

An overview of the cyclic energy balances for the DTM and CTM simulation 
runs is shown in table 5.7. 
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SIMULATED AND EXPERIMENTAL BULK VELOCITY PROFILES 

Probe at 1676 mm Reynolds no: 11706 Valensi no: 79.6 
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Figure 5.43 Discontinuous triggering mode 
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Figure 5.44 Continuous triggering mode 



Table 5.7 Comparison case cyclic energy balances. 


Triggering mode 

Net Heat 

Net Heat 

Net Kinetic 

Indicated 

Energy 

Input 

Output 

Energy Con- 

Work Done 

Balance 


m 

(j) 

version (J) 

GD 

Error (%) 

discontinuous 

0.2482 

0.7023 

0.2847 

-0.1693 

0.01016 

continuous 

0.223 

0.4347 

0.0578 

-0.1538 

0.015098 


The energy balance errors reported were realised after three simulated cycles. 
Of significance is that the net kinetic energy conversion for the DTM 
simulation is four times greater than that of the CTM simulation. This is 
attributable to the turbulence magnitude being larger in the DTM simulation 
(see figures 5.55 and 5.56). 

The mean velocity profiles at the test section are compared in figures 
5.45 to 5.47. The fluctuations apparent in figures 5.43 and 5.44 are also 
evident in figure 5.45 and are manifest through the boundary layer. The 
magnitude of the 12.5% discrepancy in the simulated and experimental maximum 
mean velocities is mainly attributable to the fluctuations in the experimental 
mean velocities. Although figures 5.46 and 5.47 are similar, the DTM 
simulation appears to correspond more closely to the experimental data over 
the first quarter of the cycle than the CTM simulation. This is confirmed by 
figures 5.48 and 5.49, particularly below 50°. As noted above, because the 
difference in core velocities is largely due tcx systematic experimental 
effects which are not included in the simulation, the impact of turbulence 
model errors on the mean velocity profiles cannot be readily determined. 

Hence, ignoring the magnitude discrepancies, figure 5.50 shows that the CTM 
simulation yields velocity profiles that parallel their experimental 
counterparts fairly well in the core flow region. However, the simulated 
velocity profiles in the boundary layer do not agree well with the 
experimental data beyond 116.5°. Returning to figures 5.48 and 5.49, the DTM 
simulation seems to yield a closer correspondence to the experimental data in 
the boundary layer than the CTM simulation. Comparing figures 5.48 to 5.50 
with figure 5.35 (showing data produced without a turbulence model), it is 
evident that inclusion of a turbulence model does yield significantly better 
correspondence between the simulated and experimental data in the core flow 
region. 

The discrepancies between the experimental and simulated boundary layer 
mean velocity profiles most likely are attributable to the turbulence 
triggering mechanism used in the simulation. The experimental axial RMS 
velocity fluctuations are shown in figure 5.51. These fluctuations are 
converted into turbulent kinetic energy (k) in figure 5.52. Strictly, this k 
cannot be compared with that simulated because the simulated k includes both 
axial and radial components while the measured k consists of an axial 
component only. The k- u> and k- c models assume that k is isotropic (that is, 
the axial and radial components are equal) which may not be true in 
oscillating flows. However, without experimental radial velocity fluctuation 
data, no further determination on this issue can be made here. 
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FIGURE 5.45 EXPERIMENTAL NEGATIVE HALF CYCLE DATA (Re=11706; Va=79.6) 




FIGURE 5.46 SIMULATION NEGATIVE HALF CYCLE DATA (Re=11706; Va=79.6) 
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FIGURE 5.47 SIMULATION NEGATIVE HALF CYCLE DATA (Re = 11706; Va=79.6) 
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Figure 5.48 Discontinuous triggering mode 
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Figure 5.49 Continuous triggering mode 
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Figure 5.50 Continuous triggering mode 




X AXIS: ANGLE [deg] 

MIN: 0.000 MAX: 180.000 
Y AXIS: RADIAL POS [mm] 
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FIGURE 5.51 EXPERIMENTAL NEGATIVE HALF CYCLE DATA (Re = 11706; Va=79.6) 
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FIGURE 5.52 EXPERIMENTAL NEGATIVE HALF CYCLE DATA (Re=11706; Va-79.6) 



With this reservation, the experimental and simulated turbulent kinetic 
energies are compared in figures 5.52 to 5.54. The effects of the 
experimental mean velocity fluctuations on the turbulent kinetic energy is 
readily apparent ih figure 5.52, which thus does not display the same 
smoothness as figures 5.53 and 5.54. The overall shape of the DTM simulation 
k surface corresponds better to the experimental data than that of the CTM 
simulation (figure 5.54). However, the DTM simulation turbulence triggering 
occurs too early compared with the experimental observations, that is, the 
discontinuous triggering level (based on bulk flow considerations) of Ng e — 
2000 is too low. It may thus be inferred that discontinuous triggering is 
physically more accurate but that a bulk flow parameter such as the Reynolds 
number is inadequate as a trigger parameter. As suggested in section 
5.6.3(c) , use of a microscopically based trigger parameter may prove to be 
physically more realistic. 

The CTM simulation depicts the existence of relatively large k over the 
entire half cycle (figure 5.54) and thus does not correspond to the physical 
behavior depicted in figure 5.52. In contrast, the DTM simulation yields a 
maximum ic 3.8 times greater than that measured while the CTM simulation and 
measured k are within 20% of each other. Assuming that the isotropic 
turbulence assumption is valid for the experimental flow, then a net maximum 
experimental k of .811 (m/s) 2 falls one-third of the way between the CTM and 
DTM simulated values. This warrants a suspicion that the CTM simulation is 
more accurate in predicting the evolution of turbulence once it has been 
triggered while the DTM simulation better reflects the triggering mechanism 
itself. Lastly, the differences in turbulence triggering revealed in figures 
5.52 to 5.54 can be seen to account for the discrepancy in the mean velocity 
profiles evident in figures 5.48 to 5.50. '* 

The full -cycle DTM and CTM simulations produce radically different 
surfaces as shown in figures 5.55 and 5.56. As there is no experimental data 
against which the simulated full-cycle k data may be directly compared, 
qualitative recourse again must be made to figure 5.38, which is the closest 
data point in Seume's thesis to the comparison case test point. Based on a 
comparison of figures 5.38, 5.55 and 5.56, it is apparent that, topologically, 
the DTM simulation produces far better correspondence with the experimental 
data than the CTM simulation. Examination of figure 5.55 reveals that the 
magnitude of the negative half-cycle turbulence (0° to 180°) is slightly less 
than the positive half-cycle turbulence, thus mimicking the trend of figure 
5.38. However, the simulated relative positive and negative half-cycle 
kinetic energies are probably in error. Figure 5.56 bears little resemblance 
to figure 5.38 revealing inverted behavior, that is, negative half-cycle k 
greater than the positive half-cycle k. Based on this comparison, it must be 
concluded that continuous triggering does not appear to be physically viable, 
at least for the two -parameter , k-u model evaluated. 
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FIGURE 5.53 SIMULATION NEGATIVE HALF CYCLE DATA (Re-11706; Va 79.6) 
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FIGURE 5.54 SIMULATION NEGATIVE HALF CYCLE DATA (Re=11706; Va=79.6) 




FIGURE 5.55 SIMULATION FULL CYCLE DATA (Re=11706; Va=79.6) 




FIGURE 5.56 SIMULATION FULL CYCLE DATA (Re=11706; Va=79.6) 



5.9 CONCLUSION 


In the context of the limited experimental data available, the results 
show that a two-parameter, k-u foldback function turbulence model 
incorporating a discontinuous triggering mode is capable of producing 
qualitatively reasonable turbulence predictions. In terms of the mean flow 
parameter predictions, which are important to Stirling machine design and 
analysis, two-parameter turbulence models of the k-u type may be adequate for 
most purposes, so making use of a more sophisticated (although rigorously 
preferable) Reynolds stress turbulence model unnecessary. 

Viewing the simulation and testing of the METR as a whole, it seems that 
the advective turbulence triggering hypothesis of Seume is the predominant 
triggering mechanism, while the local acceleration has more influence on 
controlling the evolution of the turbulence once triggering has occurred. 

This triggering mechanism is fundamentally discontinuous , requiring the 
development of a turbulence triggering parameter which is fully realizable in 
terms of the turbulence model variables. The bulk flow Reynolds number does 
not appear to be suitable because it is macroscopically rather than 
microscopically defined. The experimental and simulated data both suggest 
that a microscopically based triggering parameter is required for physical 
accuracy . 

The turbulence model development and testing reported should be viewed 
only as an exploratory first step that has highlighted the turbulence 
triggering issue as being cardinal to the successful simulation of oscillating 
turbulent flows. More experimental data and simulation development is 
necessary before the k-u turbulence model presented is suitable for 
application to the oscillating flows found in Stirling cycle machines. 
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CHAPTER 6 


CONCLUSION 


Three of the issues that have emerged during the course of this project 
as being of significance to Stirling machine simulation are: 

1. Two-dimensional component simulation of Stirling machines does not 
necessarily yield any improvement over one -dimensional simulation 
in overall performance prediction. 

2. The ability to model information propagation effects at low Mach 
numbers may be a necessary attribute for Stirling machine 
simulation codes. This determination may be made in terms of the 
characteristic number (JV C ^) which appears to be a useful 
additional non-dimensional parameter for describing Stirling 
machine oscillating flows. 

3. Turbulence triggering is a key element in implementing a 
successful turbulence model for oscillating flows. 

While simulating the METR turbulent flow field, the need for reliable 
oscillating turbulent flow experimental data spanning a broad range of 
boundary conditions has become apparent. As advective boundary conditions 
have a major, if not overriding, significance on the turbulence triggering 
mechanism in oscillating flows, particular care should be taken to quantify 
these boundary conditions completely. However, it must be recognised that 
turbulence models validated using experimental data generated with non- 
Stirling advective boundary conditions may not be immutably applicable to 
Stirling, machines themselves because of the inherently chaotic nature of the 
process. Thus, ultimately, such turbulence models may have to be 
validated against turbulent flow data gathered in Stirling machines directly. 

Although two-dimensional component simulation is clearly of value in 
understanding the flow details of Stirling machines, from the work carried out 
here, it is not yet evident that such component simulation provides any better 
prediction of overall machine performance than one - dimensional simulation. 

From a design perspective, two-dimensional component simulation thus may serve 
best as a means of generating improved turbulence correlations for use in 
existing one - dimensional codes. Nevertheless, the continuing evolution of 
computer hardware will eventually make end-to-end, two-dimensional simulation 
of Stirling machines a cost-effective reality. 

Use of qualifying benchmark tests for Stirling machine simulation codes 
may prove to be an effective means of ensuring that these codes meet a minimum 
standard of validity, particularly if they are used to design novel (and 
expensive) Stirling hardware. Iberall's (Ib50) experimentally validated 
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(Wa 65 , G068) analytic solution of the transmission line problem is an ideal 
candidate for a qualifying benchmark test. 

As a general proposition, numerical methods are an important facet of 
Stirling machine simulation to the extent that the numerical method itself 
emulates a physical process. Thus, simulation errors may be a result of the 
physical implications of a particular numerical method, rather than just its 
numerical accuracy implications. In this light, an isolated investigation of 
pressure- linking as a legitimate numerical algorithm for simulating low 
characteristic number compressible flows is warranted. 
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